Regression

Das Vorgehen in diesem Notebook orientiert sich am Data Science Lifecycle. Es wird eine Modellpräzision von > 70% gewünscht. Das bedeutet, dass der R-squared Wert des besten Modells >0,7 sein soll.

Python Setup

Zunächst werden mit der Import Methode wichtige Python Bibliotheken importiert.

import pandas as pd
import numpy as np 

import seaborn as sns
sns.set_theme(style='ticks')

import matplotlib.pyplot as plt
%matplotlib inline

import statsmodels.formula.api as smf
import statsmodels.api as sm
import statsmodels.gam.api as smg
from statsmodels.compat import lzip
from statsmodels.stats.outliers_influence import variance_inflation_factor
from statsmodels.tools.tools import add_constant
from patsy import dmatrix


from sklearn.compose import ColumnTransformer
from sklearn.compose import make_column_selector as selector
from sklearn.pipeline import Pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler, OneHotEncoder
from sklearn.model_selection import train_test_split
import sklearn.metrics as skm
from sklearn.linear_model import LinearRegression
from sklearn.linear_model import Lasso
from sklearn.preprocessing import SplineTransformer
from sklearn.linear_model import Ridge
from sklearn.pipeline import make_pipeline

import plotly.express as px
import plotly.graph_objects as go
import plotly.offline as pyo
# Set notebook mode to work in offline
pyo.init_notebook_mode()

Datenaufbereitung

Daten laden

Importieren der Daten aus der CSV in Github und Erstellung eines Pandas Dataframes zur Bearbeitung und Verwendung der Daten in diesem Notebook.

GIT = "https://raw.githubusercontent.com/jan-kirenz/project-annikas428/main/"
DATA = "project_data.csv"
TOKEN = "?token=GHSAT0AAAAAABPURMBOEW2WS3JEBXGI5TPAYPJRBIQ"

df_initial = pd.read_csv(GIT + DATA + TOKEN)

Daten inspizieren

Im ersten Schritt ist es wichtig, sich einen Überblick über die Daten zu verschaffen. Begonnen wird mit einer Anzeige einiger Datensätze.

df_initial
longitude latitude housing_median_age total_rooms total_bedrooms population households median_income median_house_value ocean_proximity price_category
0 -122.23 37.88 41.0years 880 129.0 322 126 8.3252 452600.0$ NEAR BAY above
1 -122.22 37.86 21.0 7099 1106.0 2401 1138 8.3014 358500.0 NEAR BAY above
2 -122.24 37.85 52.0 1467 190.0 496 177 7.2574 352100.0 NEAR BAY above
3 -122.25 37.85 52.0 1274 235.0 558 219 5.6431 341300.0 NEAR BAY above
4 -122.25 37.85 52.0 1627 280.0 565 259 3.8462 342200.0 NEAR BAY above
... ... ... ... ... ... ... ... ... ... ... ...
20635 -121.09 39.48 25.0 1665 374.0 845 330 1.5603 78100.0 INLAND above
20636 -121.21 39.49 18.0 697 150.0 356 114 2.5568 77100.0 INLAND above
20637 -121.22 39.43 17.0 2254 485.0 1007 433 1.7000 92300.0 INLAND above
20638 -121.32 39.43 18.0 1860 409.0 741 349 1.8672 84700.0 INLAND above
20639 -121.24 39.37 16.0 2785 616.0 1387 530 2.3886 89400.0 INLAND above

20640 rows × 11 columns

Beim Sichten der Datensätze fällt auf, dass im ersten Datensatz zwei Unstimmigkeiten im Vergleich zu den anderen Datensätzen vorliegen:
median_house_value ist nur hier mit der Währung Dollar angegeben und hinter housing_median_age steht der Zusatz “years”.

Nun folgt Pandas Dataframe Info Methode. Diese zeigt die Variablen und deren Datentypen, sowie die Anzahl der non-null values je Spalte an.

df_initial.info()
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 20640 entries, 0 to 20639
Data columns (total 11 columns):
 #   Column              Non-Null Count  Dtype  
---  ------              --------------  -----  
 0   longitude           20640 non-null  float64
 1   latitude            20640 non-null  float64
 2   housing_median_age  20640 non-null  object 
 3   total_rooms         20640 non-null  int64  
 4   total_bedrooms      20433 non-null  float64
 5   population          20640 non-null  int64  
 6   households          20640 non-null  int64  
 7   median_income       20640 non-null  float64
 8   median_house_value  20640 non-null  object 
 9   ocean_proximity     20640 non-null  object 
 10  price_category      20640 non-null  object 
dtypes: float64(4), int64(3), object(4)
memory usage: 1.7+ MB

Die Variablen housing_median_age und und median_house_value wurden als Datentyp object definiert. Dies muss genauer untersucht werden, da bei einer Jahreszahl bzw. einem monetären Wert ein numerischer Datentyp (int oder float) angemessener ist. Die Ursache hierfür ist auf die Zusatzangaben “$” und “years” zurückzuführen.
Die Variablen ocean_proximity und price_category sind keine numerischen sondern kategoriale Datentypen. Diese lassen sich jeweils in einer internen Rangordnung anordnen. Sie sind daher ordinal. Der richtige Datentyp hierfür ist category. Die Transformation der Datentypen erfolgt in Kapitel “Datentransformation”.
Neben den Datentypen wird auch die Anzahl der non-null values je Variable angezeigt. Hier haben alle Variablen eine Anzahl von 20640 non-null values. Lediglich die Variable total_bedrooms liegt unter diesem Wert. Dies lässt darauf schließen, dass hier einige Datensätze keinen Wert haben. Dies gilt es im Folgenden genauer zu untersuchen.

# show missing values -if present - in yellow
sns.heatmap(df_initial.isnull(), 
            yticklabels=False,
            cbar=False, 
            cmap='viridis');

print(df_initial.isnull().sum())
longitude               0
latitude                0
housing_median_age      0
total_rooms             0
total_bedrooms        207
population              0
households              0
median_income           0
median_house_value      0
ocean_proximity         0
price_category          0
dtype: int64
_images/regression_13_1.png

Wie bereits vermutet, enthält die Variable total_bedrooms 207 null values. Die Abbildung visualisiert diese nicht befüllten Datensätze und deren Verteilung im Datenset in gelb. Die null-values sind zufällig im Datenset verteilt. Da das Datenset über 20000 Datensätze enthält, werden die Zeilen mit den null-values im Kapitel Transform Data entfernt. Alternativ wäre auch eine Belegung der null-values beispielsweise mit dem Durchschnittswert der Variable möglich. Darauf wird hier jedoch verzichtet, genug andere Datensätze vorliegen.

Datentransformation

Entfernung der null values

# drop rows with missing values
df = df_initial.dropna()
print(df.isnull().sum())
longitude             0
latitude              0
housing_median_age    0
total_rooms           0
total_bedrooms        0
population            0
households            0
median_income         0
median_house_value    0
ocean_proximity       0
price_category        0
dtype: int64

Entfernung unwichtiger Features

Das Feature price_category wird entfernt, da es sich direkt aus der Variablen median_house_value ableitet und daher nicht zur Berechnung des Wertes genutzt werden darf. Dieses Feature wird im Notebook Classification als abhängige Variable genutzt.

#drop irrelevant features
df = df.drop(['price_category'], axis=1)

Konvertierung von Datentypen

Die mithilfe der info() Methode identifizierten Features, welche eine Datentyp Konvertierung benötigen, werden nun bearbeitet.

# Convert data types Ocean_Proximity and Price_Category to Data Type "Category" as they are categorical variables (ordinal)
df.ocean_proximity = df.ocean_proximity.astype('category')

# Housing_Median_Age and Median_House_Value are actually int64 values, but pandas did not get this because of signs like '€'
df.housing_median_age = df.housing_median_age.str.replace('years', '', regex=False)
df.median_house_value = df.median_house_value.str.replace('$', '', regex=False)
df.housing_median_age = pd.to_numeric(df.housing_median_age, errors='coerce').astype('int64')
df.median_house_value = pd.to_numeric(df.median_house_value, errors='coerce').astype('float64')
df.total_bedrooms = df.total_bedrooms.astype('int64')
df.info()
df.head()
<class 'pandas.core.frame.DataFrame'>
Int64Index: 20433 entries, 0 to 20639
Data columns (total 10 columns):
 #   Column              Non-Null Count  Dtype   
---  ------              --------------  -----   
 0   longitude           20433 non-null  float64 
 1   latitude            20433 non-null  float64 
 2   housing_median_age  20433 non-null  int64   
 3   total_rooms         20433 non-null  int64   
 4   total_bedrooms      20433 non-null  int64   
 5   population          20433 non-null  int64   
 6   households          20433 non-null  int64   
 7   median_income       20433 non-null  float64 
 8   median_house_value  20433 non-null  float64 
 9   ocean_proximity     20433 non-null  category
dtypes: category(1), float64(4), int64(5)
memory usage: 1.6 MB
longitude latitude housing_median_age total_rooms total_bedrooms population households median_income median_house_value ocean_proximity
0 -122.23 37.88 41 880 129 322 126 8.3252 452600.0 NEAR BAY
1 -122.22 37.86 21 7099 1106 2401 1138 8.3014 358500.0 NEAR BAY
2 -122.24 37.85 52 1467 190 496 177 7.2574 352100.0 NEAR BAY
3 -122.25 37.85 52 1274 235 558 219 5.6431 341300.0 NEAR BAY
4 -122.25 37.85 52 1627 280 565 259 3.8462 342200.0 NEAR BAY

Mithilfe der string replace Methode werden die zuvor identifizierten Zeichen “years” und “$” entfernt. Anschließend ist eine Transformation der Features in integer bzw. float Werte möglich. Für housing_median_age wird integer gewählt, da die Jahreszahlen als Ganze Zahlen also discrete Werte angegeben sind.
Für median_house_value wird der Datentyp float benötigt, da die Ergebnisse einer linearen Regression stehts continous sind.

Hinzufügen neuer Features

df = df.assign(people_per_household=lambda df: df.population / df.households)
df = df.assign(rooms_per_household=lambda df: df.total_rooms / df.households)
df = df.assign(bedrooms_per_people=lambda df: df.total_bedrooms / df.population)
df = df.assign(bedrooms_per_rooms=lambda df: df.total_bedrooms / df.total_rooms)

Diese Features werden dem Datenset hinzugefügt, in der Hoffnung, dass diese für die Verwendung in den Modellen nützlich sein könnten. Untersucht wird dies im Kapitel EDA.

Datensplit

Split des Datensets in 80% Trainingsdaten und 20% Testdaten. Die Testdaten werden anschließend bis zur finalen Evaluation der besten Modelle nicht verwendet.

train_dataset = df.sample(frac=0.8, random_state=0) 
test_dataset = df.drop(train_dataset.index)

train_dataset
longitude latitude housing_median_age total_rooms total_bedrooms population households median_income median_house_value ocean_proximity people_per_household rooms_per_household bedrooms_per_people bedrooms_per_rooms
14185 -117.08 32.70 37 2176 418 1301 375 2.8750 98900.0 NEAR OCEAN 3.469333 5.802667 0.321291 0.192096
6125 -117.91 34.11 20 3158 684 2396 713 3.5250 153000.0 <1H OCEAN 3.360449 4.429173 0.285476 0.216593
14095 -117.10 32.75 11 2393 726 1905 711 1.3448 91300.0 NEAR OCEAN 2.679325 3.365682 0.381102 0.303385
14359 -117.22 32.74 52 1260 202 555 209 7.2758 345200.0 NEAR OCEAN 2.655502 6.028708 0.363964 0.160317
18004 -121.99 37.29 32 2930 481 1336 481 6.4631 344100.0 <1H OCEAN 2.777547 6.091476 0.360030 0.164164
... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
6168 -117.93 34.09 35 1891 353 1093 382 4.0167 165500.0 <1H OCEAN 2.861257 4.950262 0.322964 0.186674
4174 -118.21 34.10 36 2000 533 1234 535 3.7437 241700.0 <1H OCEAN 2.306542 3.738318 0.431929 0.266500
3593 -118.48 34.24 32 2621 412 1285 414 6.6537 267600.0 <1H OCEAN 3.103865 6.330918 0.320623 0.157192
5366 -118.39 34.04 44 1873 286 635 283 5.5951 461300.0 <1H OCEAN 2.243816 6.618375 0.450394 0.152696
17879 -122.00 37.40 17 5121 1017 1470 968 2.9706 81300.0 <1H OCEAN 1.518595 5.290289 0.691837 0.198594

16346 rows × 14 columns

Deskriptive Statistik

Kennzahlenanalyse numerischer Variablen

Mithilfe der describe() Methode werden die wichtigsten Kennzahlen numerischer Variablen (bspw. Durchschnitt, min und max Werte) ausgegeben.

#summary statistics for all numerical columns
train_dataset.describe().T
count mean std min 25% 50% 75% max
longitude 16346.0 -119.564154 2.002618 -124.350000 -121.790000 -118.490000 -118.000000 -114.470000
latitude 16346.0 35.630318 2.138574 32.550000 33.930000 34.250000 37.710000 41.950000
housing_median_age 16346.0 28.664505 12.556764 1.000000 18.000000 29.000000 37.000000 52.000000
total_rooms 16346.0 2622.235776 2169.548287 11.000000 1448.000000 2119.000000 3120.750000 39320.000000
total_bedrooms 16346.0 535.281659 418.469078 3.000000 296.000000 432.500000 644.000000 6445.000000
population 16346.0 1416.087055 1103.842065 3.000000 784.250000 1164.000000 1711.000000 28566.000000
households 16346.0 496.758167 379.109535 3.000000 280.000000 408.000000 600.000000 6082.000000
median_income 16346.0 3.869337 1.902228 0.499900 2.555675 3.533200 4.744225 15.000100
median_house_value 16346.0 206916.154411 115676.394484 14999.000000 119300.000000 179700.000000 265900.000000 500001.000000
people_per_household 16346.0 3.089529 11.525259 0.692308 2.428850 2.816514 3.282584 1243.333333
rooms_per_household 16346.0 5.423747 2.160962 0.846154 4.443697 5.230149 6.051703 62.422222
bedrooms_per_people 16346.0 0.403370 0.223233 0.000670 0.315480 0.372646 0.443446 8.750000
bedrooms_per_rooms 16346.0 0.213188 0.057891 0.100000 0.175591 0.203263 0.239668 1.000000

Der Durchschnitt (mean) der Variablen housing_median_age über alle Datensätze beträgt 28,64 Jahre. Die Standardabweichung, also die durchschnittliche Abweichung der Werte vom Durchschnitt, beträgt 12,56 Jahre. Das Minimum liegt bei 1 Jahr, das Maximum bei 52 Jahren. 25% der Datensätze liegen unter 18 Jahren (1. Quartil). 50% unter 29 Jahren (2. Quartil) und 75% unter 37 Jahren (3. Quartil).
Im Durchschnitt beträgt der Wert people_per_household 3,09. Die durchschnittliche Abweichung von diesem Wert beträgt 11,53. Das Minimum liegt bei 0,69 und das Maximum bei 1243 Personen pro Haushalt. 25% der Datensätze liegen zwischen 0,69 und 2,43. 50% unter 2,82 und 75% der Datensätze unter 3,28 Personen pro Haushalt.

print(train_dataset.quantile(0.99, interpolation='nearest'))
longitude                 -116.290000
latitude                    40.630000
housing_median_age          52.000000
total_rooms              11021.000000
total_bedrooms            2208.000000
population                5723.000000
households                1955.000000
median_income               10.557500
median_house_value      500001.000000
people_per_household         5.286920
rooms_per_household         10.500000
bedrooms_per_people          0.941606
bedrooms_per_rooms           0.407472
Name: 0.99, dtype: float64

Die Quantile() Methode gibt für jede Variable einen Wert an, unter welchem (in diesem Beispiel) 99% der Datensätze liegen.
Beispielsweise liegen 99% der Datensätze unter einer Anzahl total_rooms in Höhe von 11021.
99% der Distrikte haben eine population unter 5723 Einwohnern.
Im Nachfolgenden werden einige mögliche Ausreißer identifiziert.

train_dataset.loc[train_dataset.people_per_household > 70]
longitude latitude housing_median_age total_rooms total_bedrooms population households median_income median_house_value ocean_proximity people_per_household rooms_per_household bedrooms_per_people bedrooms_per_rooms
9172 -118.59 34.47 5 538 98 8733 105 4.2391 154600.0 INLAND 83.171429 5.123810 0.011222 0.182156
19006 -121.98 38.32 45 19 5 7460 6 10.2264 137500.0 INLAND 1243.333333 3.166667 0.000670 0.263158
16669 -120.70 35.32 46 118 17 6532 13 4.2639 350000.0 NEAR OCEAN 502.461538 9.076923 0.002603 0.144068
3364 -120.51 40.41 36 36 8 4198 7 5.5179 67500.0 INLAND 599.714286 5.142857 0.001906 0.222222

Mehr als 70 Leute in einem Haushalten wirken sehr unrealistisch. 99% der Datensätze liegen unter einem Wert von 5,29 Personen pro Haushalt. Daher werden diese Datensätze als Ausreißer identifiziert.

train_dataset.loc[train_dataset.bedrooms_per_people > 8]
longitude latitude housing_median_age total_rooms total_bedrooms population households median_income median_house_value ocean_proximity people_per_household rooms_per_household bedrooms_per_people bedrooms_per_rooms
11862 -121.25 40.27 25 958 245 28 16 2.625 67500.0 INLAND 1.75 59.875 8.75 0.255741

Auch dieser Wert liegt weit über dem 99% Perzentil der Varialben bedrooms_per_people und wird daher als Ausreiser entfernt.

train_dataset.loc[train_dataset.rooms_per_household > 55]
longitude latitude housing_median_age total_rooms total_bedrooms population households median_income median_house_value ocean_proximity people_per_household rooms_per_household bedrooms_per_people bedrooms_per_rooms
12447 -114.49 33.97 17 2809 635 83 45 1.6154 87500.0 INLAND 1.844444 62.422222 7.650602 0.226059
11862 -121.25 40.27 25 958 245 28 16 2.6250 67500.0 INLAND 1.750000 59.875000 8.750000 0.255741
1913 -120.06 39.01 19 2967 528 112 48 4.0714 437500.0 INLAND 2.333333 61.812500 4.714286 0.177958

Auch diese Datensätze liegen weit über dem Wert des 99% Perzentils und werden daher aus dem Datenset entfernt.

#drop people per household outliers
train_dataset = train_dataset.drop([3364, 9172, 16669, 19006])
#drop bedroooms per people outliers
train_dataset = train_dataset.drop([11862])
#drop rooms per household_outliers
train_dataset = train_dataset.drop([1913, 12447])
#summary statistics for all numerical columns
train_dataset.describe().T
count mean std min 25% 50% 75% max
longitude 16339.0 -119.564115 2.002470 -124.350000 -121.790000 -118.490000 -118.000000 -114.470000
latitude 16339.0 35.629562 2.138070 32.550000 33.930000 34.250000 37.710000 41.950000
housing_median_age 16339.0 28.664973 12.555983 1.000000 18.000000 29.000000 37.000000 52.000000
total_rooms 16339.0 2622.903544 2169.632035 11.000000 1449.000000 2120.000000 3121.500000 39320.000000
total_bedrooms 16339.0 535.416978 418.477283 3.000000 296.000000 433.000000 644.000000 6445.000000
population 16339.0 1415.032315 1100.485866 3.000000 785.000000 1164.000000 1711.000000 28566.000000
households 16339.0 496.956301 379.069221 3.000000 280.000000 408.000000 600.500000 6082.000000
median_income 16339.0 3.869002 1.901829 0.499900 2.555600 3.532900 4.744150 15.000100
median_house_value 16339.0 206925.109248 115665.635455 14999.000000 119300.000000 179700.000000 265900.000000 500001.000000
people_per_household 16339.0 2.941847 1.128881 0.692308 2.428863 2.816514 3.282400 63.750000
rooms_per_household 16339.0 5.413424 2.023751 0.846154 4.443662 5.230061 6.051152 52.848214
bedrooms_per_people 16339.0 0.402250 0.202972 0.014902 0.315552 0.372651 0.443430 7.925926
bedrooms_per_rooms 16339.0 0.213190 0.057897 0.100000 0.175587 0.203262 0.239667 1.000000

Analyse kategorialer Variablen

# summary statistics for all categorical columns
train_dataset.describe(include=['category']).transpose()
count unique top freq
ocean_proximity 16339 5 <1H OCEAN 7216

Die Variable ocean_proximity besitzt 5 Ausprägungen, wovon die Ausprägung <1H OCEAN mit einer Anzhl von 7216 am häufigsten vorkommt.

train_dataset['ocean_proximity'].value_counts()
<1H OCEAN     7216
INLAND        5215
NEAR OCEAN    2107
NEAR BAY      1796
ISLAND           5
Name: ocean_proximity, dtype: int64

Hier wird die Häufigkeit einzelner Ausprägungen von ocean_proximity dargestellt. Die Ausprägung ISLAND kommt über das gesamte Datenset mit über 20000 Datensätzen nur fünfmal vor. In Relation zur Gesamtanzahl der Datensätzen ist dies sehr wenig, wodurch die Aussagekraft dieser Ausprägung vermutlich eher gering ist. Daher werden die Datensätze mit Ausprägung ISLAND mit der Ausprägung NEAR OCEAN vereint.

train_dataset.loc[train_dataset.ocean_proximity == 'ISLAND']
longitude latitude housing_median_age total_rooms total_bedrooms population households median_income median_house_value ocean_proximity people_per_household rooms_per_household bedrooms_per_people bedrooms_per_rooms
8318 -118.48 33.43 29 716 214 422 173 2.6042 287500.0 ISLAND 2.439306 4.138728 0.507109 0.298883
8316 -118.32 33.33 52 2127 512 733 288 3.3906 300000.0 ISLAND 2.545139 7.385417 0.698499 0.240715
8315 -118.33 33.34 52 2359 591 1100 431 2.8333 414700.0 ISLAND 2.552204 5.473318 0.537273 0.250530
8314 -118.32 33.35 27 1675 521 744 331 2.1579 450000.0 ISLAND 2.247734 5.060423 0.700269 0.311045
8317 -118.32 33.34 52 996 264 341 160 2.7361 450000.0 ISLAND 2.131250 6.225000 0.774194 0.265060
# write out a dict with the mapping of old to new
remap_cat_dict = {
    'NEAR OCEAN': 'NEAR OCEAN',
    'NEAR BAY': 'NEAR BAY',
    'ISLAND': 'NEAR OCEAN',
    'INLAND': 'INLAND',
    '<1H OCEAN': '<1H OCEAN' }

train_dataset.ocean_proximity = train_dataset.ocean_proximity.map(remap_cat_dict).astype('category')

Exploratory Data Anaysis

Dieses Kapitel dient der ausführlichen Datenanalyse der Testdaten. Hierfür werden vor allem geeignete Visualisierungen in Form von Diagrammen verwendet.

erster Überblick

fig = go.Figure(data=go.Scattergeo(
        lon = train_dataset['longitude'],
        lat = train_dataset['latitude'],
        text = 'median house value in district: ' + train_dataset['median_house_value'].astype(str) + '$',
        marker_color = train_dataset['median_house_value'], marker_size=2, marker_colorscale = 'Inferno', marker_colorbar= {'title': 'Median House Value'}
        ))

fig.update_layout(
        title = 'Median House Value per District in California',
        geo_scope='usa',
    )

pyo.iplot(fig)

Die Datensätze stammen alle aus Kalifornien, USA.
Die Werte Longitude und Latitude unterscheiden sich daher kaum, wie sich auch bereits im Kapitel deskriptive Statistik gezeigt hat. Daher sind sie als Feature eher ungeeignet.
Allerdings scheint die Entfernung zum Ozean eine wichtige Rolle im Bezug auf den median_house_value zu spielen. ocean_proximity scheint daher als Feature interessant und sollte genauer untersucht werden.

train_dataset = train_dataset.drop(['longitude', 'latitude'], axis = 1)
sns.pairplot(train_dataset[['housing_median_age', 'total_rooms', 'total_bedrooms', 'population', 'households', 'median_house_value']])
<seaborn.axisgrid.PairGrid at 0x209e49bc850>
Error in callback <function flush_figures at 0x00000209E2DE5040> (for post_execute):

Die hier abgebildeten, unabhängigen Variablen weisen eine starke Korrelation untereinander auf (bis auf housing_median_age).
In Bezug auf die abhängige Variable median_house_value zeigen die Plots zunächst keinen starken Zusammenhang.
Interessant sieht das Histogramm von median_house_value aus, da es am äußeren rechten Rand nochmal hohe Häufigkeiten haben. Dieses Histogramm wird später nochmal genauer analysisiert.

sns.pairplot(train_dataset[['median_income', 'people_per_household', 'rooms_per_household', 'bedrooms_per_people', 'bedrooms_per_rooms', 'median_house_value']]);
_images/regression_56_0.png

Zwischen median_income und median_house_value ist ein linearer Zusammenhang erkennbar. Auch das Histogramm von median_income zeigt eine hohe Säule am rechten Rand.
Auch für die selbst definierten Variablen ist es schwierig, ein linearen Zusammenhang zu median_house_value zu erkennen.

fig = sns.displot(data= train_dataset, x='median_house_value')
_images/regression_58_0.png

Die Verteilung von Median_House_value ist nicht normalverteilt, rechts im hohen Preissegment ist die Häufigkeit nochmal extrem hoch. Über 500000 Dollar exisiteren keine weiteren Hauswerte. Möglicherweise wurde der Datensatz mit Werten über 500000 Dollar abgeschnitten und für alle höheren Werte, der Wert 500000 gesetzt.

fig = sns.displot(data= train_dataset, x='median_income')
_images/regression_60_0.png

Womöglich wurde auch die Variable median_income nach oben abgeschnitten.

Korrelationsanalyse

# Create correlation matrix for numerical variables
corr_matrix = train_dataset.corr()
corr_matrix['median_house_value'].sort_values(ascending=False)
median_house_value      1.000000
median_income           0.692099
rooms_per_household     0.176503
total_rooms             0.130828
housing_median_age      0.105534
bedrooms_per_people     0.077487
households              0.063096
total_bedrooms          0.047606
population             -0.026634
people_per_household   -0.178408
bedrooms_per_rooms     -0.253025
Name: median_house_value, dtype: float64
# Make a pretty heatmap

# Use a mask to plot only part of a matrix
mask = np.zeros_like(corr_matrix)
mask[np.triu_indices_from(mask)]= True

# Change size
plt.subplots(figsize=(11, 15))

# Build heatmap with additional options
heatmap = sns.heatmap(corr_matrix, 
                      mask = mask, 
                      square = True, 
                      linewidths = .5,
                      cmap = 'coolwarm',
                      cbar_kws = {'shrink': .6,
                                'ticks' : [-1, -.5, 0, 0.5, 1]},
                      vmin = -1,
                      vmax = 1,
                      annot = True,
                      annot_kws = {"size": 10})
_images/regression_64_0.png

Laut der Korrelationsmatrix einigen sich vor allem folgende Variablen:
median_income weist mit einem Wert von 0,69 die höchste Korrelation auf, danach folgen bedrooms_per_rooms mit einer negativen Korrelation in Höhe von -0,25 und anschließend people_per_household ebenfalls mit einer negativen Korrelation in Höhe von -0,18.
Danach folgt rooms_per_household, allerdings kann diese Variable nicht in Kombination mit bedrooms_per_rooms oder people_per_household genutzt werden, da sie sich aus total_rooms und households zusammensetzt. Diese Variablen werden bereits in bedrooms_per_rooms und people_per_household verwendet, daher ist die zusätzliche Verwendung von rooms_per_household nicht sinnvoll. Die Korrelation der Variablen untereinander ist in diesem Fall zu hoch mit einem Wert von -0,51 bzw. 0,6.
housing_median_age könnte mit einer Korrelation von 0,11 noch einen positiven Einfluss auf das Modell haben.

Analyse der kategorialen Variable Ocean_Proximity

sns.boxenplot(data=train_dataset, x='ocean_proximity', y='median_house_value')
<AxesSubplot:xlabel='ocean_proximity', ylabel='median_house_value'>
_images/regression_67_1.png
sns.stripplot(data=train_dataset, x='ocean_proximity', y='median_house_value', size=2)
<AxesSubplot:xlabel='ocean_proximity', ylabel='median_house_value'>
_images/regression_68_1.png

Die Ausprägungen NEAR BAY und NEAR OCEAN unterscheiden sich hinsichtlich median_house_value kaum. Daher können diese beiden Ausprägungen zusammengefasst werden.

sns.displot(data=train_dataset, x='median_house_value', hue='ocean_proximity', kind='kde')
<seaborn.axisgrid.FacetGrid at 0x230d8270f40>
_images/regression_70_1.png

Auch in der Häufigkeitsverteilung zeigen NEAR BAY und NEAR OCEAN kaum Unterschiede.

# write out a dict with the mapping of old to new
remap_cat_dict = {
    'NEAR OCEAN': 'NEAR COAST',
    'NEAR BAY': 'NEAR COAST',
    'INLAND': 'INLAND',
    '<1H OCEAN': '<1H OCEAN' }

train_dataset.ocean_proximity = train_dataset.ocean_proximity.map(remap_cat_dict).astype('category')
sns.boxplot(data=train_dataset, x='ocean_proximity', y='median_house_value')
<AxesSubplot:xlabel='ocean_proximity', ylabel='median_house_value'>
_images/regression_73_1.png

Auch NEAR COAST und <1H OCEAN zeigen eine ähnliche Verteilung des Medians und der Quartile. Das 3. Quartil von NEAR COAST zieht sich jedoch noch um ca. 50000$ weiter nach oben, als das von <1H OCEAN, daher bleibt die Differenzierung dieser beiden Ausprägungen noch bestehen.
Bei der Ausprägung INLAND sind viele Ausreiser nach oben zu erkennen.

Analyse numerischer Variablen

# check relationship with a joint plot
sns.jointplot(x="median_income", y="median_house_value", hue="ocean_proximity", data=train_dataset);
_images/regression_76_0.png
# check relationship with a joint plot
sns.lmplot(x="median_income", y="median_house_value", data=train_dataset);
_images/regression_77_0.png

Wie bereits aus der Korrelationsmatrix zu schließen, ist ein linearer Zusammenhang zwischen median_income und median_house_value erkennbar. Die “Linie” mit Datenpunkten bei 500000 Dollar ergibt sich aus der anzunehmenden Begrenzung der median_house_value_ auf 500000 Dollar.
Auch innerhalb der einzelenden Ausprägungen von ocean_proximity ist eine lineare Korrelation von median_income und median_house_value erkennbar.

sns.relplot(data=train_dataset, x='total_rooms', y='median_house_value', hue='ocean_proximity', col='ocean_proximity')
<seaborn.axisgrid.FacetGrid at 0x193868f47c0>
_images/regression_79_1.png

Ein linearer Zusammenhang von total_rooms zu median_house_value ist nicht wirklich zu erkennen. Die Variable ist daher für eine lineare Regression eher ungeeignet.

# Plot with Plotly Express
fig = px.scatter(x=train_dataset['bedrooms_per_rooms'], y=train_dataset['median_house_value'], opacity=0.65, 
                trendline='ols', trendline_color_override='darkred')
pyo.iplot(fig)
fig = px.scatter(x=train_dataset['people_per_household'], y=train_dataset['median_house_value'], opacity=0.65, 
                trendline='ols', trendline_color_override='darkred')
pyo.iplot(fig)

bedrooms_per_rooms und people_per_household zeigen einen negativen linearen Zusammenhang. Bereits in der Korrelationsmatrix war der negative Zusammenhang dieser beiden Variablen zu median_house_value sichtbar. Die Variablen werden nachfolgend in die linearen Regressionsmodelle aufgenommen.

Feature Auswahl

Auf Basis der in der deskriptiven Statistik und EDA gewonnen Erkenntnisse, sollen nun geeignete Features für die Modellierung ausgewählt werden.
Diese Auswahl wird auf Basis der sogenannten Filter Methode getroffen.
Diese eine von drei Möglichkeiten bei der Feature Selektion. Dabei werden beispielsweise auf Basis einer Korrelationsmatrix die besten Feature ausgewählt und nur diese dem Modell zur Verfügung gestelltn.
Folgende Features wurden als geeignet identifiziert:

  1. median_income

  2. ocean_proximity

  3. bedrooms_per_rooms

  4. people_per_household

  5. housing_median_age (eventuell)

Modellierung

Statsmodel

Model 1 - OLS Regression

Als erstes Modell wird statsmodels OLS Regression verwendet.
Im ersten Modell werden lediglich die Feature 1-4 verwendet.

# Fit Model
lm = smf.ols(formula='median_house_value ~ median_income + ocean_proximity + bedrooms_per_rooms + people_per_household' , data=train_dataset).fit()
# Full summary
lm.summary()
OLS Regression Results
Dep. Variable: median_house_value R-squared: 0.623
Model: OLS Adj. R-squared: 0.623
Method: Least Squares F-statistic: 5407.
Date: Wed, 19 Jan 2022 Prob (F-statistic): 0.00
Time: 12:34:59 Log-Likelihood: -2.0569e+05
No. Observations: 16339 AIC: 4.114e+05
Df Residuals: 16333 BIC: 4.114e+05
Df Model: 5
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept 3.394e+04 4569.346 7.428 0.000 2.5e+04 4.29e+04
ocean_proximity[T.INLAND] -6.785e+04 1408.255 -48.179 0.000 -7.06e+04 -6.51e+04
ocean_proximity[T.NEAR COAST] 1.557e+04 1423.620 10.934 0.000 1.28e+04 1.84e+04
median_income 4.323e+04 404.027 107.007 0.000 4.24e+04 4.4e+04
bedrooms_per_rooms 3.009e+05 1.3e+04 23.197 0.000 2.75e+05 3.26e+05
people_per_household -1.376e+04 496.667 -27.713 0.000 -1.47e+04 -1.28e+04
Omnibus: 5010.041 Durbin-Watson: 1.972
Prob(Omnibus): 0.000 Jarque-Bera (JB): 27166.062
Skew: 1.373 Prob(JB): 0.00
Kurtosis: 8.689 Cond. No. 128.


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
lm.params
Intercept                         33939.853081
ocean_proximity[T.INLAND]        -67848.787400
ocean_proximity[T.NEAR COAST]     15566.550456
median_income                     43233.585628
bedrooms_per_rooms               300850.020547
people_per_household             -13763.891665
dtype: float64

Die summary() Methode berechnet die wichtigsten Kennzahlen des Modells.
R-squared ist eine Kennzahl zur Beurteilung der Anpassungsgüte eines Modells und nimmt einen Wert zwischen 0% und 100% an. Als Bezugsbasis wird der Durchschnitt verwendet.
lm erreicht einen R-squared Wert in Höhe von 62,3%. Somit können 62,3% der median_house_value Streuung
In statistischen Methoden wird meist lieber der adjusted R-squared genutzt. Da dieser die degrees of freedom miteinbezieht, lässt sich durch adj. R-squared einen besseren Rückschluss auf die Gesamtpopulation der Datensätze ziehen. Dieser ist meistens etwas schlechter als R-squared. In diesem Modell sind beide Werte gleich hoch.
Mithilfe der Kennzahlen AIC und BIC soll ein overfitting des Modells verhindert werden. Für jedes zusätzliche Feature wird eine “Bestrafung” in die Kennzahlen miteinberechnet. Je kleiner die Werte, desto besser. In diesem Modell sind die Werte jedoch relativ hoch.
F-statistics ist die Menge der systematischen Varianz (MSM) geteilt durch die Menge der unsystematischen Varianz (MSR). Die Kennzahl gibt an, in welcher Höhe das Modell die Ergebnisausgabe der abhängigen Variable verbessert hat, verglichen zu Ungenauigkeit im Modell. Je höher F-statistics, desto besser. Ein F-statistics Wert in Höhe von 5407, wie in Modell lm, ist schon sehr gut.
Der Y-Achsenschnittpunkt des Modells liegt bei 33939,85. Die Steigung der einzelnen Features ist in der Spalte coef abzulesen. Zu interpretieren sind diese folgendermaßen: Vergrößert sich der Wert in Variable median_income um 1 und alle anderen Werte bleiben konstant, so ist das Ergebnis um 43233,59 größer als zuvor.
Die p-values der Features liegen alle bei 0,00. Die Nullhypothese (ß=0 => Feature hat keinen Einfluss auf die abhängige Variable) kann demnach für alle Features verworfen werden.
Die p-values für Jarque-Bera und Omnibus liegen ebenfalls bei 0,00. Die Nullhypothese, die Fehler seien normalverteilt, wird in diesem Falle verworfen. Es wird also angenommen, die Residuals seien nicht normalverteilt. Bei einer Datensetgröße > 50, kann dies aufgrund des Central Limit Theorems jedoch vernachlässigt werden.
Der Durbin-Watson Wert liegt bei 1,97. Somit kann angenommen werden, dass die Fehler des Modells keine Korrelation aufweisen.

Regression Diagnostics
fig = sm.graphics.influence_plot(lm, criterion="cooks")
fig.tight_layout(pad=1.0)
_images/regression_94_0.png

Mithilfe des Influence Plots werden mögliche Ausreißer als blaue Punkte in der Größe ihres Einflusses identifiziert. Dies können Werte mit hoher Abweichung in ihrem Y-Wert (outliers) oder X-Wert (high leverage) sein. In einer linearen Regression können solche Ausreise einen großen Einfluss auf das Gesamtmodell haben und somit die Performance des gesamten Modells erheblich verschlechtern. Als Faustregel werden Datenpunkte mit einer Cooks Distance in Höhe von 4/n (wobei n die Anzahl der Datensätze ist) als mögliche Ausreißer identifiziert.

# obtain Cook's distance 
lm_cooksd = lm.get_influence().cooks_distance[0]

# get length of df to obtain n
n = len(train_dataset["median_income"])
print('n = ', n)

# calculate critical d
critical_d = 4/n
print('Critical Cooks distance:', critical_d)

# identification of potential outliers with leverage
out_d = lm_cooksd > critical_d

# output potential outliers with leverage
print(train_dataset.index[out_d], "\n", 
    lm_cooksd[out_d])
n =  16339
Critical Cooks distance: 0.0002448130240528796
Int64Index([15778,  6776, 11511,  7263, 15723, 14756, 15687,  4638, 14806,
            15815,
            ...
            15809,  8866,   199, 18498,    59, 10749,  8857, 13154,  5470,
            16642],
           dtype='int64', length=934) 
 [1.49490398e-03 2.68414884e-04 4.78133170e-04 5.02665948e-04
 5.94624792e-04 1.45661981e-02 6.06459084e-04 9.89584645e-04
 7.26161610e-04 8.30878122e-04 1.36260242e-03 9.64095843e-04
 4.86783348e-04 3.28695576e-04 3.75058682e-04 7.93946113e-04
 5.82170067e-04 2.59400713e-04 5.38430141e-04 4.34597064e-04
 2.58773707e-04 4.39852205e-04 3.24124357e-04 3.81088866e-04
 2.95276469e-04 8.44927027e-04 4.83475406e-02 5.21391591e-04
 4.36852748e-04 4.06798688e-04 3.74228212e-04 2.75462275e-04
 8.37885379e-04 4.91036338e-04 1.65550557e-03 2.28062263e-03
 1.11362268e-03 5.40836931e-04 5.06680630e-04 2.80028526e-03
 3.84825684e-04 4.88676566e-04 8.68935319e-04 2.51788726e-03
 6.76056994e-04 3.19987726e-04 3.94500298e-04 8.87693727e-04
 3.88788562e-04 4.33792026e-04 3.43876382e-04 4.28843898e-04
 7.54230456e-04 3.50021703e-04 4.80369640e-04 5.30319765e-04
 4.59051144e-04 6.00454845e-04 3.80966666e-04 4.31199706e-04
 4.03345609e-04 8.57174195e-04 4.11096007e-04 4.63221951e-04
 1.23566980e-03 2.50497036e-04 4.82794226e-04 5.26525139e-04
 5.69167729e-04 3.28935485e-04 6.97236942e-04 1.22926028e-03
 7.37881318e-04 2.93996391e-04 3.92057182e-04 4.70614674e-04
 1.17905418e-03 5.23060291e-04 4.86889004e-04 5.40614399e-04
 1.95440918e-03 2.90055671e-04 9.44897532e-04 7.33532196e-04
 6.47125775e-04 3.18823752e-04 1.06578506e-03 3.58388250e-04
 8.73670560e-04 1.10228431e-03 7.93461030e-04 2.61856821e-04
 2.98663391e-04 5.60713494e-04 8.96529290e-04 2.49653913e-04
 1.07462907e-03 3.66163173e-04 3.27982973e-04 6.74318828e-04
 1.07600427e-03 3.12835081e-04 9.41739908e-04 1.06035239e-03
 3.24403292e-04 2.32717050e-03 3.33725982e-04 3.31713778e-04
 2.45370019e-04 5.19796814e-04 2.52998567e-04 3.15019545e-04
 7.17730666e-04 4.95060839e-04 5.84604429e-04 6.76035191e-04
 2.62979595e-04 2.78405274e-04 2.64614306e-04 3.43991908e-04
 8.28547805e-04 3.96502354e-04 5.72061896e-04 9.03304907e-04
 9.88962274e-04 3.88578207e-04 6.03733226e-04 4.42201426e-04
 5.84262789e-04 3.02906624e-04 6.67400127e-04 4.21938365e-04
 3.18350166e-02 2.50731877e-04 7.56862933e-04 3.26525900e-04
 2.91886314e-04 4.94062027e-04 2.62415163e-04 7.40021481e-04
 7.74038819e-04 6.64113946e-04 6.00653167e-04 2.59384025e-04
 5.68021977e-04 5.08871869e-04 3.45793984e-04 3.16681659e-04
 3.59376409e-04 2.75726762e-04 5.71869366e-04 2.55386427e-04
 2.78580587e-04 4.80446565e-04 6.29032892e-04 2.68844084e-04
 3.15777697e-04 4.18502209e-04 2.66244188e-04 9.72580016e-04
 3.59936346e-04 1.26716166e-03 4.61208595e-03 2.72426608e-04
 3.31337333e-04 6.13153738e-04 3.94343301e-04 3.35424565e-04
 5.11719688e-04 2.87763886e-04 2.48788511e-04 5.29875909e-04
 2.52462604e-04 6.18585440e-04 2.93302623e-04 2.59391598e-04
 2.57900725e-04 2.67878340e-04 7.54047576e-04 7.90100598e-04
 3.51351563e-04 8.86252032e-04 1.18258562e-03 2.64701347e-04
 3.06120700e-03 5.21490540e-04 3.94491160e-04 2.69647319e-04
 2.65030920e-04 1.24489084e-03 6.75547522e-04 3.70755646e-03
 9.64789182e-04 2.46399968e-04 3.96264268e-04 2.99146468e-04
 5.65433172e-04 2.98636628e-03 3.31239647e-04 3.41387971e-04
 2.69398818e-04 6.40829474e-04 3.15785273e-03 2.61734223e-04
 3.60648977e-04 2.75598565e-04 4.26681256e-04 1.74170377e-02
 3.74987588e-03 2.79744243e-03 8.47589712e-04 1.34331315e-03
 4.85411722e-04 4.38351271e-04 3.92230730e-04 2.69865195e-04
 6.13534479e-04 6.00850811e-04 2.56881788e-04 3.00847558e-04
 2.60229616e-04 3.71622289e-04 4.01045487e-04 4.40508578e-04
 9.71587095e-04 3.40572615e-04 2.68812469e-04 3.70078254e-04
 6.10179623e-04 3.02152403e-04 2.67996019e-04 7.44260318e-04
 3.13479879e-04 8.19067610e-02 5.12803460e-04 3.93272586e-04
 2.19823559e-03 3.53653351e-04 4.34626944e-04 1.39726947e-03
 1.18617182e-03 5.01736586e-04 5.30094203e-04 3.33634176e-04
 4.18932862e-02 9.42177142e-04 4.18773467e-04 6.03680488e-04
 3.71553393e-04 3.74286350e-04 1.23669485e-03 4.65135400e-04
 2.69268542e-04 4.13374391e-04 3.53588165e-04 9.81277198e-04
 2.54137328e-04 1.35995340e-03 3.11246270e-03 4.31574146e-04
 7.86358637e-04 2.69481369e-04 2.98644303e-04 5.35187123e-04
 2.51539827e-04 3.13073701e-04 9.94563742e-04 3.79299795e-04
 2.72704600e-04 2.61286323e-04 2.94840250e-04 2.61668973e-03
 2.65733675e-04 3.11871682e-04 1.16617792e-03 3.22950323e-03
 3.59684839e-04 1.61885963e-03 4.58416247e-04 2.48030609e-04
 3.42383137e-04 4.48068631e-04 6.04910702e-04 8.04956576e-04
 3.11900750e-04 2.50644499e-04 3.76788326e-04 5.23889548e-04
 3.98319955e-04 2.47144472e-04 3.04257523e-03 2.97617992e-04
 2.70824454e-04 5.76938158e-04 7.24114626e-04 3.09400255e-03
 5.36686516e-03 4.42145065e-04 3.02082029e-04 2.63000798e-04
 3.88936728e-04 3.38298305e-04 4.44264842e-04 2.71776042e-04
 2.56975665e-04 4.40378588e-04 9.46502316e-04 4.44336623e-04
 3.68241357e-03 8.19569279e-04 7.25573623e-03 4.63097728e-04
 3.86542144e-03 1.79765644e-03 3.13796286e-04 1.85849209e-03
 3.63517011e-04 3.05689630e-04 4.46005083e-04 1.86968485e-03
 4.21898380e-04 6.28339395e-04 2.74402091e-04 5.88073246e-04
 2.32709199e-03 4.52837607e-04 1.50790793e-03 5.81240516e-04
 1.64359329e-03 5.80674834e-03 6.55261768e-04 1.73632609e-03
 3.32072802e-04 4.10811940e-04 1.04001291e-03 5.98785387e-04
 2.79515437e-04 6.23027502e-04 4.23275012e-04 1.38618438e-03
 2.71617878e-04 2.39674723e-03 3.51234700e-04 3.86146581e-04
 4.07595499e-04 4.20280178e-03 3.69799077e-04 2.58852189e-04
 1.00625914e-03 2.14059605e-03 5.24085589e-04 5.86927033e-04
 2.76194574e-04 6.77176762e-04 6.22637307e-04 6.22956497e-04
 1.12480619e-03 4.16146778e-04 3.70163512e-04 7.00008802e-04
 3.59599089e-04 3.82851777e-04 4.55744907e-04 1.10934093e-03
 4.52817097e-04 2.87570210e-04 2.52329760e-04 7.02649092e-04
 2.59392154e-04 4.46540428e-04 5.57525212e-04 4.63626144e-04
 3.65869117e-04 5.78289495e-03 6.15022887e-04 3.35508687e-04
 2.65061093e-04 2.02798819e-03 2.49284001e-04 3.23955147e-04
 8.10369090e-04 3.80530786e-04 4.15664656e-04 2.95094039e-04
 8.47132464e-04 2.85414761e-04 2.84178883e-04 2.65532294e-04
 2.71743672e-04 3.20261144e-04 4.67532383e-04 3.33304384e-04
 2.52355395e-04 1.28385469e-03 2.87227641e-04 7.19836520e-04
 1.24407306e-03 3.60495943e-04 2.99981365e-04 4.65103270e-04
 5.46220633e-04 4.25385196e-04 1.69057212e-03 9.95180258e-04
 3.27190396e-04 2.85787247e-04 3.28116249e-04 7.50507763e-04
 4.31914472e-04 3.33314939e-04 3.96174751e-04 4.33313938e-04
 8.35842338e-04 3.91664204e-04 5.53686911e-04 4.22118493e-04
 2.49932882e-03 7.81329389e-04 2.75321278e-03 2.71498687e-04
 3.92285936e-04 1.11281331e-03 2.80607976e-04 4.96095023e-03
 4.21027426e-04 3.23368164e-04 4.14159128e-04 3.77467629e-03
 6.07651838e-04 7.18666140e-04 4.53756366e-04 3.41057171e-03
 2.77265109e-04 2.55933137e-03 3.99882274e-04 5.34894447e-04
 5.20700768e-04 2.52589779e-04 8.76764590e-04 6.01984749e-04
 1.01172240e-03 4.34738639e-04 4.86857764e-03 4.44271558e-04
 3.77835060e-04 4.17273122e-04 2.49841045e-04 8.57922366e-04
 1.01138406e-03 3.17770892e-04 2.79425670e-04 2.87924766e-04
 2.75321438e-04 2.16673451e-03 5.56504395e-04 3.85884793e-03
 6.33212441e-04 6.84666057e-04 3.21755993e-04 5.23263229e-04
 4.07589874e-04 1.36587639e-03 8.41514470e-04 6.10850040e-04
 2.87242598e-04 3.28588051e-01 3.71667189e-04 2.79196718e-03
 4.09579459e-04 4.37807917e-04 8.51212481e-04 4.86433645e-04
 9.30857717e-04 5.01044984e-04 2.47930811e-04 2.93683726e-03
 3.59552309e-04 2.57442888e-04 1.45592560e-03 3.45294621e-04
 7.20625186e-04 2.50894360e-04 7.58404371e-04 3.63506552e-04
 2.32150175e-03 3.57160051e-04 2.82166230e-04 4.05308960e-04
 1.05660916e-03 2.36117472e-03 2.55300114e-04 2.99426708e-04
 2.48755421e-04 3.47614085e-04 5.35389815e-04 2.28333284e-03
 2.69689693e-04 1.83993869e-02 2.78134600e-04 2.94646617e-04
 6.32580902e-04 3.37268213e-04 8.49437828e-04 4.70592108e-04
 2.65122577e-04 7.87353472e-04 2.65151373e-04 3.40147104e-04
 1.47373921e-03 3.58812729e-04 6.71821996e-03 3.17543351e-04
 3.61536099e-04 9.04605585e-03 5.89081995e-04 7.59534268e-04
 2.64907542e-04 2.85244797e-03 6.43833067e-03 3.42696744e-04
 3.53282313e-04 3.02325255e-04 8.96809026e-04 2.57985742e-04
 8.29591390e-04 2.56536173e-03 3.98997462e-04 3.80381504e-04
 4.04846515e-04 5.53919005e-04 1.72091124e-03 3.41707784e-04
 4.54876276e-04 8.50177014e-04 3.12956299e-04 1.32323097e-03
 3.46349028e-04 3.41588528e-04 5.65716565e-04 5.51966438e-04
 4.49538242e-04 5.01573516e-04 2.67487102e-04 7.86724676e-04
 2.66813565e-04 2.81069561e-04 3.72920665e-04 2.67001014e-04
 2.45976115e-04 8.53312926e-04 3.47462338e-04 3.12518224e-04
 3.37634912e-04 6.71193185e-03 2.80483340e-04 1.05535803e-03
 5.05727691e-04 5.49198981e-04 3.92816712e-04 3.26462785e-03
 3.78658807e-04 1.07133515e-03 2.58061845e-04 6.28308148e-04
 5.99432217e-04 4.08464699e-04 7.76299893e-03 3.32998562e-04
 7.75864814e-04 2.43445894e-03 1.29339473e-03 3.88447003e-04
 2.72155160e-04 1.16799456e-03 6.99554601e-04 3.41654550e-04
 2.60497617e-04 6.14377217e-04 1.12340141e-03 9.52380819e-04
 3.28633882e-04 2.66993087e-03 7.63434605e-04 2.72831557e-04
 4.40328575e-04 3.27503237e-04 9.76400663e-04 8.92877629e-04
 8.15770615e-04 3.34734715e-04 6.30237071e-04 3.76652904e-04
 3.53213140e-04 4.80060218e-04 3.10904734e-04 5.28275910e-04
 2.94585293e-04 2.14279533e-03 3.85897937e-04 5.67212170e-04
 2.67990920e-04 4.60165209e-04 7.28971845e-04 8.93621528e-04
 7.40300668e-04 5.67891373e-04 2.98114027e-03 5.06945094e-04
 2.90940363e-04 7.08720513e-04 2.58135066e-04 2.93029122e-04
 2.81941867e-04 9.58097289e-04 5.82149611e-04 4.05186998e-04
 3.70767144e-03 2.78477079e-03 2.89663179e-04 1.20228109e-03
 4.54930522e-04 2.70933886e-04 4.04780848e-04 7.11395439e-04
 3.32404421e-04 2.79054479e-04 3.51763013e-03 3.66814342e-04
 6.17607727e-04 1.92440482e-03 7.40987676e-04 6.02952970e-04
 2.90803016e-04 2.87937670e-04 4.13298793e-04 6.02532557e-04
 7.03667841e-04 2.85359319e-04 4.43822131e-04 4.48611300e-04
 1.84253469e-03 2.45013962e-04 5.37653017e-04 3.52961018e-04
 6.78910376e-04 1.02483048e-03 1.04860911e-03 1.49712177e-03
 2.69398997e-04 4.40503340e-04 2.55594132e-04 1.19645327e-03
 2.53989647e-04 4.47037845e-04 3.51513607e-04 4.08201711e-04
 1.50074578e-03 8.30344714e+00 5.65625002e-04 2.46168492e-04
 5.33248865e-04 1.22758173e-03 1.01208124e-03 4.25612245e-04
 6.64096985e-04 2.47406741e-04 2.70997639e-04 2.97073012e-04
 3.13625041e-04 1.40881899e-03 3.02927343e-04 7.98328272e-04
 4.70190253e-04 3.24999226e-04 3.06698993e-04 3.93282175e-04
 4.77023978e-04 6.30043037e-04 9.15404010e-04 3.07118690e-04
 3.75683525e-04 2.69553212e-04 6.54946966e-04 2.95261234e-04
 3.85318203e-03 5.86215552e-04 3.24299295e-04 9.44380999e-04
 4.08059729e-04 1.78452639e+00 2.78842121e-03 2.93019660e-04
 5.81763831e-04 2.46515742e-04 2.67421853e-04 2.46015947e-04
 2.59353083e-04 3.30689996e-04 3.10118580e-04 5.66627749e-04
 3.76468283e-04 5.23602207e-04 3.17035827e-04 4.32554132e-04
 2.88423909e-03 3.18177800e-04 2.47097467e-04 8.97296353e-04
 4.11799006e-03 3.76872367e-04 5.61415054e-04 4.44857116e-04
 2.99348646e-03 2.51524523e-04 5.19464921e-04 5.07202394e-04
 4.18908548e-03 4.46062016e-04 1.23169781e-03 4.41564831e-04
 4.31238292e-04 3.78250481e-04 2.59072590e-04 2.57755138e-04
 2.56018898e-04 5.76159655e-04 1.67865668e-03 2.98192987e-04
 2.79505548e-04 5.43658186e-04 2.55414133e-04 2.57925230e-04
 7.50001815e-04 2.64633710e-04 3.71728892e-04 5.23219628e-04
 3.08927700e-04 2.78169248e-04 3.03316115e-04 3.29103640e-04
 3.50804495e-04 2.47784628e-04 3.43024543e-04 8.04772883e-04
 3.68602838e-04 3.63458034e-03 3.61478735e-03 4.11608446e-04
 3.75786430e-04 2.62410102e-04 4.64802687e-04 2.63496724e-04
 5.33857567e-04 2.77274533e-04 5.01078667e-04 3.42090208e-04
 7.95906573e-04 3.48992531e-04 2.74019023e-04 3.82462266e-04
 3.49499382e-04 3.60329524e-03 4.98543051e-04 4.58754170e-04
 2.99287784e-04 2.90066791e-04 5.11005725e-04 8.67868112e-04
 3.03929160e-04 2.85321238e-04 1.03605731e-02 6.19344809e-04
 2.46206722e-04 2.61928905e-03 2.89126514e-04 4.04317282e-04
 3.14789439e-04 5.41032937e-04 4.68125716e-04 2.82483142e-04
 3.63130780e-04 7.50205131e-04 2.66739351e-04 4.18751814e-04
 7.22603412e-04 5.51362034e-04 7.91641879e-04 4.36816848e-04
 3.79457580e-04 4.46525265e-04 3.12461430e-03 2.76876967e-04
 2.58483722e-04 3.53791532e-03 2.97204974e-04 2.62999419e-04
 3.97663321e-04 5.85841186e-03 9.09442576e-04 2.83333635e-04
 5.26047782e-04 9.12294548e-04 3.20830967e-04 3.99222484e-04
 6.95444999e-04 1.67166712e-03 3.23984602e-04 4.27941932e-04
 4.11748278e-04 1.31224780e-03 1.07865097e-02 2.73367535e-04
 3.07708117e-04 3.23280459e-04 2.69735178e-04 3.60040135e-04
 2.85055681e-03 2.86618779e-04 2.49277549e-04 4.60375511e-04
 2.96411485e-04 6.52882093e-04 4.90617969e-04 3.95001225e-04
 3.21764867e-04 2.92752399e-04 4.63077330e-04 5.00309549e-04
 3.32037995e-04 3.20749075e-04 2.54231674e-04 4.00198976e-04
 2.54940788e-04 9.70628458e-04 2.94543782e-04 3.36519756e-04
 2.93415300e-04 2.49376829e-04 6.95317814e-04 3.61124579e-04
 4.69950482e-04 2.91806508e-04 2.85707017e-03 2.38613707e-03
 2.63286803e-04 2.88384114e-04 1.85190590e-03 9.93274311e-04
 2.94307710e-04 2.91103242e-04 3.76488541e-04 7.12226233e-04
 4.58230854e-04 3.20857506e-03 3.27931387e-04 6.53579046e-04
 2.80913497e-03 4.47200022e-04 4.98029532e-04 2.59734059e-04
 5.67699033e-04 2.97946489e-04 2.51796775e-04 7.55961243e-04
 5.70450094e-04 2.80284751e-04 7.34322115e-04 3.82171751e-04
 3.13375643e-04 6.55482693e-04 2.78098379e-04 2.96805132e-04
 6.86610201e-04 3.28259515e-04 4.71620905e-04 1.04241011e-03
 6.19521175e-04 2.76714681e-04 3.06755289e-04 4.62806729e-04
 3.69218926e-04 5.44672688e-03 3.70524290e-04 3.07698834e-03
 8.89940345e-04 2.84599116e-04 6.65044184e-04 3.12831265e-04
 3.56531287e-04 3.72030288e-04 3.27770091e-04 3.54258528e-04
 2.57750234e-04 3.87354535e-04 4.36528660e-03 2.96099098e-04
 9.38024385e-01 4.32715339e-04 2.37162395e-03 4.54278443e-04
 5.15313448e-04 1.24983495e-03 3.78456415e-04 1.24498717e-03
 5.83661630e-04 4.62907874e-04 4.10149561e-04 3.86241023e-04
 2.83869049e-04 3.01381616e-03 1.83772898e-03 3.77798003e-04
 3.53370743e-04 3.07065505e-04 2.98605919e-04 2.73322810e-04
 8.09264921e-04 3.14676919e-04 4.63606496e-04 4.30723878e-04
 5.55919618e-04 5.29171324e-04 4.66093611e-04 2.10386753e-03
 3.33256599e-04 5.02692880e-04 2.87112698e-04 3.74359189e-04
 3.12368308e-04 8.57238075e-04 2.80246652e-04 3.21797874e-04
 1.01113123e-03 5.10846722e-03]

In diesem Modell werden 934 Datensätze mit einer Cooks Distance > 4/n als mögliche Ausreißer identifiziert. Beispielsweise könnte versucht werden, den Einfluss dieser Datensätze auf das Modell durch Datentransformation (bspw. durch Verwendung des Logarithmus) zu verringern. Allerdings hilft diese Methode der Modellperformance genau so oft, wie sie ihr schadet. Daher wird hier auf eine Datentransformation verzichtet und die Datenpunkte stattdessen aus dem Trainingsset entfernt.

#drop outliers identified by Cook's distance
train_dataset2 = train_dataset.drop(train_dataset.index[out_d])
train_dataset2.info()
<class 'pandas.core.frame.DataFrame'>
Int64Index: 15405 entries, 14185 to 17879
Data columns (total 12 columns):
 #   Column                Non-Null Count  Dtype   
---  ------                --------------  -----   
 0   housing_median_age    15405 non-null  int64   
 1   total_rooms           15405 non-null  int64   
 2   total_bedrooms        15405 non-null  int64   
 3   population            15405 non-null  int64   
 4   households            15405 non-null  int64   
 5   median_income         15405 non-null  float64 
 6   median_house_value    15405 non-null  float64 
 7   ocean_proximity       15405 non-null  category
 8   people_per_household  15405 non-null  float64 
 9   rooms_per_household   15405 non-null  float64 
 10  bedrooms_per_people   15405 non-null  float64 
 11  bedrooms_per_rooms    15405 non-null  float64 
dtypes: category(1), float64(6), int64(5)
memory usage: 1.4 MB
fig = sns.scatterplot(data=train_dataset2, x='median_income', y='median_house_value')
_images/regression_99_0.png

Im Vergleich zum ursprünglichen Trainingsset, welches im Kapitel EDA ausführlich analysiert wurde, fehlen hier einige X- und Y-Ausreißer. Dadurch hat sich die waagrechte Punktelinie bei 500000$ medin_house_value verkleinert.

fig = sns.displot(data=train_dataset2, x='median_house_value')
_images/regression_101_0.png

Durch die Entfernung der Ausreißer hat sich das Histogramm von median_house_value etwas weiter normalisiert.

# Fit Model without outliers
lm = smf.ols(formula='median_house_value ~ median_income + ocean_proximity + bedrooms_per_rooms + people_per_household', data=train_dataset2).fit()
lm.summary()
OLS Regression Results
Dep. Variable: median_house_value R-squared: 0.741
Model: OLS Adj. R-squared: 0.741
Method: Least Squares F-statistic: 8832.
Date: Wed, 19 Jan 2022 Prob (F-statistic): 0.00
Time: 12:37:35 Log-Likelihood: -1.8925e+05
No. Observations: 15405 AIC: 3.785e+05
Df Residuals: 15399 BIC: 3.786e+05
Df Model: 5
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept 3.596e+04 4371.548 8.226 0.000 2.74e+04 4.45e+04
ocean_proximity[T.INLAND] -6.127e+04 1092.791 -56.071 0.000 -6.34e+04 -5.91e+04
ocean_proximity[T.NEAR COAST] 6831.6716 1113.689 6.134 0.000 4648.711 9014.633
median_income 4.676e+04 362.691 128.920 0.000 4.6e+04 4.75e+04
bedrooms_per_rooms 3.536e+05 1.2e+04 29.400 0.000 3.3e+05 3.77e+05
people_per_household -2.532e+04 576.811 -43.895 0.000 -2.64e+04 -2.42e+04
Omnibus: 901.834 Durbin-Watson: 1.985
Prob(Omnibus): 0.000 Jarque-Bera (JB): 1114.776
Skew: 0.584 Prob(JB): 8.50e-243
Kurtosis: 3.612 Cond. No. 153.


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

Durch die Entfernung der Ausreißer haben sich R-squared, adjusted R-squared und F-statistics deutlich verbessert.
AIC und BIC sind auch etwas kleiner geworden und haben sich somit verbessert.
Die p-values der Features bleiben unverändert.
Auch an den p-values für Jarque-Bera und Omnibus hat sich kaum etwas verändert.

fig = sm.graphics.plot_partregress_grid(lm)
fig.tight_layout(pad=0.2)
_images/regression_105_0.png

Die Partial Regression Plots visualisieren die Relation einer unabhängigen Variablen in Verbindung mit den anderen unabhängigen Variablen im Modell in Bezug auf die abhängige Variable.
Die Plots zeigen eine negative Abhängigkeit von T.INLAND und peope_per_household auf median_house_value. Median_income und bedrooms_per_rooms scheinen eine positive Abhängigkeit auf median_house_value zu haben. Gemäß dem Diagramm hat die Variable ocean_proximity[T.NEAR COAST] keinen Einfluss auf median_income. Dies steht im Gegensatz zum p-value 0,00 dieser Variablen, nach dessen Interpretation die Variable einen Einfluss auf die abhängige Variable hat.

fig = sm.graphics.plot_regress_exog(lm, "ocean_proximity[T.NEAR COAST]")
fig.tight_layout(pad=1.0)
_images/regression_107_0.png

Auch diese Plots lassen darauf schließen, dass die Variable _ocean_proximity[T.NEAR COAST] das Modell nicht verbessert. Sowohl Partial Regression Plot als auch CCPR Plot visualisieren eine Nicht-Linearität zwischen der Variablen und median_house_value.
Die Residuals scheinen dagegen gleichmäßig verteilt zu sein.

fig = sm.graphics.plot_regress_exog(lm, "ocean_proximity[T.INLAND]")
fig.tight_layout(pad=1.0)
_images/regression_109_0.png

Die negative Korrelation ist sowohl im Partial Regression Plot als auch im CCPR Plot zu erkennen. Die Residuals sind gleichmäßig verteilt und lassen daher auf Homoskedastizität rückschließen. Das Diagramm zu Y and Fitted vs. X ist aufgrund der vielen Datenpunkte schwer zu analysieren.

fig = sm.graphics.plot_regress_exog(lm, "people_per_household")
fig.tight_layout(pad=1.0)
_images/regression_111_0.png

Auch hier ist eine negative Abhängigkeit der abhängigen zur unabhängigen Variablen zu erkennen. Linearität scheint also vorzuliegen. Allerdings ist im Residual Plot keine konstante Varianz der Fehlerstreuung zu erkennen. Mit steigendem X-Wert verringert sich die Fehlerstreuung und liegt nur noch im positiven Bereich. Mit zunehmender people_per_household wird tendenziell ein niedrigerer median_house_value mit dem Modell hervorgesagt, als tatsächlich vorliegt. Dies spiegelt sich auch im Plot zu Y and Fitted vs. X wieder. Bei dieser Variablen könnte Heteroskedastizität vorlieren. Dies gilt es mithilfe des Breusch-Pagan Lagrange Multiplier Tests genauer zu untersuchen:

name = ['Lagrange multiplier statistic', 'p-value', 'f-value', 'f p-value']
test = sm.stats.het_breuschpagan(lm.resid, lm.model.exog)
lzip(name, test)
[('Lagrange multiplier statistic', 1043.0254994601257),
 ('p-value', 2.906700494490131e-223),
 ('f-value', 223.66770899896423),
 ('f p-value', 3.329323077460106e-231)]

Da alle Werte über 0,05 liegen, kann die Nullhypothese, Homoskedastizität, akzeptiert werden. Somit liegt keine Hetereoskedastizität (Korrelation der Fehlerterms) vor.

fig = sm.graphics.plot_regress_exog(lm, "median_income")
fig.tight_layout(pad=1.0)
_images/regression_115_0.png

Für median_income ist die Linearität im Partial Regression Plot und CCPR Plot deutlich zu erkennen. Jedoch nimmt auch hier die Streuung der Residuals im höheren X-Wertebereich ab. Heteroskedastizität wurde mithilfe des Breusch-Pagan Lagrange Multiplier Tests bereits ausgeschlossen.

Für das nächste Modell werden die Ausprägungen NEAR COST und <1H OCEAN zusammengefasst:

# write out a dict with the mapping of old to new
remap_cat_dict = {
    'NEAR COAST': 'NEAR COAST',
    'INLAND': 'INLAND',
    '<1H OCEAN': 'NEAR COAST' }

train_dataset.ocean_proximity = train_dataset.ocean_proximity.map(remap_cat_dict).astype('category')

Model 2 - OLS Regression

Nachfolgend wird ein weiteres Modell mit zusätzlicher Verwendung der Variablen housing_median_age trainiert.

# Fit Model
lm2 = smf.ols(formula='median_house_value ~ median_income + ocean_proximity + bedrooms_per_rooms + people_per_household + housing_median_age', data=train_dataset).fit()
lm2.summary()
OLS Regression Results
Dep. Variable: median_house_value R-squared: 0.632
Model: OLS Adj. R-squared: 0.631
Method: Least Squares F-statistic: 5598.
Date: Wed, 19 Jan 2022 Prob (F-statistic): 0.00
Time: 12:57:18 Log-Likelihood: -2.0552e+05
No. Observations: 16339 AIC: 4.110e+05
Df Residuals: 16333 BIC: 4.111e+05
Df Model: 5
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept -5.941e+04 4251.872 -13.972 0.000 -6.77e+04 -5.11e+04
ocean_proximity[T.NEAR COAST] 6.572e+04 1340.461 49.026 0.000 6.31e+04 6.83e+04
median_income 4.432e+04 403.597 109.805 0.000 4.35e+04 4.51e+04
bedrooms_per_rooms 2.979e+05 1.28e+04 23.242 0.000 2.73e+05 3.23e+05
people_per_household -1.439e+04 487.311 -29.524 0.000 -1.53e+04 -1.34e+04
housing_median_age 1009.9369 45.939 21.984 0.000 919.891 1099.983
Omnibus: 5235.108 Durbin-Watson: 1.968
Prob(Omnibus): 0.000 Jarque-Bera (JB): 34340.227
Skew: 1.379 Prob(JB): 0.00
Kurtosis: 9.545 Cond. No. 766.


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

R-squared und adj. R-squared haben sich im Vergleich zum ersten Modell verbessert.
Auch F-statistics hat sich durch die Entfernung der Ausprägung von ocean_proximity und Aufnahme des Features median_housing_age verbessert im Vegleich zu Modell 1 ohne Entfernung der Cook’s Distance Outliers.
Alle übrigen Kennzahlen haben sich kaum verändert.

# obtain Cook's distance 
lm_cooksd = lm2.get_influence().cooks_distance[0]

# get length of df to obtain n
n = len(train_dataset["median_income"])
print('n = ', n)

# calculate critical d
critical_d = 4/n
print('Critical Cooks distance:', critical_d)

# identification of potential outliers with leverage
out_d = lm_cooksd > critical_d

# output potential outliers with leverage
print(train_dataset.index[out_d], "\n", 
    lm_cooksd[out_d])
n =  16339
Critical Cooks distance: 0.0002448130240528796
Int64Index([15778,  6776, 10005, 11511,  7263, 15723, 14756, 15687,  4638,
            14806,
            ...
            15809,  8866,   199, 18498,    59, 10749,  8857, 13154,  5470,
            16642],
           dtype='int64', length=855) 
 [1.33657972e-03 3.30853996e-04 2.57049009e-04 3.39660405e-04
 6.01129532e-04 3.88602976e-04 1.90478556e-02 6.27488571e-04
 6.74476388e-04 4.45045388e-04 8.27915922e-04 1.30029194e-03
 1.22608518e-03 3.77876614e-04 5.13746812e-04 7.71608777e-04
 5.53086958e-04 2.60583344e-04 3.05990146e-04 5.38968511e-04
 3.38650064e-04 2.54349021e-04 4.32750022e-04 3.44212356e-04
 3.61317235e-04 3.67082550e-04 2.50921062e-04 2.64010646e-04
 5.64470925e-02 4.50645263e-04 4.32630903e-04 3.73737786e-04
 8.94441250e-04 5.17859249e-04 1.67765900e-03 2.15661985e-03
 1.05909753e-03 2.68151181e-04 6.29731152e-04 5.30310328e-04
 5.08429929e-03 5.19828279e-04 9.65158028e-04 4.65438401e-03
 4.30201917e-04 3.11323002e-04 2.63976791e-04 9.17685336e-04
 2.92775090e-04 3.04159208e-04 4.29591546e-04 3.16585850e-04
 5.73111270e-04 3.40214397e-04 3.69651789e-04 4.12626701e-04
 1.19689965e-03 3.24624515e-04 6.25756905e-04 3.70903135e-04
 3.22827592e-04 3.95367319e-04 1.05559881e-03 2.65758810e-04
 2.47592532e-04 8.69636078e-04 2.51911954e-04 4.41768820e-04
 6.76831782e-04 5.65039902e-04 3.33787249e-04 7.45292365e-04
 1.02617996e-03 5.39873094e-04 4.40784778e-04 4.95448437e-04
 9.26999608e-04 2.77340040e-04 1.84706601e-03 2.91896086e-04
 8.69739979e-04 2.17140974e-03 3.09705912e-04 2.47578772e-04
 7.29639650e-04 7.24422015e-04 6.79907575e-04 3.33754238e-04
 3.15815108e-04 1.02444650e-03 3.00226007e-04 8.52987194e-04
 1.19880125e-03 9.80927157e-04 5.65469194e-04 9.45727023e-04
 1.21371291e-03 3.47389825e-04 6.50989809e-04 4.49479516e-04
 1.65423758e-03 3.89116048e-04 9.94555835e-04 7.00191096e-04
 2.18971913e-03 2.85484826e-04 2.49718013e-04 2.47484819e-04
 2.81402090e-04 5.88963038e-04 3.45832307e-04 6.95567672e-04
 6.35200588e-04 4.26719218e-04 2.96215870e-04 2.51307265e-04
 2.55813701e-04 3.63345137e-04 8.04986478e-04 2.78317071e-04
 3.96319393e-04 8.63988528e-04 1.32463248e-03 2.87769237e-04
 4.52357194e-04 3.93632745e-04 3.92028465e-04 4.05636189e-04
 3.50348707e-04 5.69254539e-04 4.26315065e-04 4.36699075e-02
 7.33412362e-04 5.49398317e-04 7.26014229e-04 1.33508397e-03
 7.01538018e-04 2.47242706e-04 5.71897176e-04 4.94021524e-04
 3.90467300e-04 3.25255208e-04 3.04205545e-04 3.20955852e-04
 7.00539166e-04 2.65241798e-04 3.52632314e-04 4.16500906e-04
 4.33311859e-04 4.41249366e-04 2.86018112e-04 9.48344105e-04
 2.47757062e-04 1.36943368e-03 5.11072467e-03 4.05563886e-04
 3.26238220e-04 2.53555665e-04 4.54670905e-04 2.61021149e-04
 3.67450043e-04 4.35513441e-04 3.45787798e-04 6.76666771e-04
 1.23889590e-03 7.43737556e-04 2.79617027e-04 2.58099362e-04
 1.38076353e-03 3.30948644e-03 3.16852303e-04 3.51331295e-04
 1.46597653e-03 7.91184137e-04 4.22485707e-03 9.02085835e-04
 3.95820079e-04 2.73525513e-04 3.67308349e-04 1.23677595e-03
 4.90126981e-03 3.65851977e-04 5.99874448e-04 2.58979835e-04
 4.28871739e-03 2.63559446e-04 4.61483732e-04 2.46606867e-02
 3.76826150e-03 3.96818944e-03 5.81099998e-04 1.47960745e-03
 3.60684593e-04 3.97478169e-04 2.64012779e-04 3.59647903e-04
 1.17153655e-03 2.51225904e-04 3.34967441e-04 3.21499259e-04
 2.71827753e-04 6.23453463e-04 1.02455371e-03 2.88236000e-04
 3.28846792e-04 4.54854733e-04 3.02734668e-04 2.56153363e-04
 1.02535591e-03 7.74844873e-02 2.95057072e-04 2.74284063e-04
 2.12112604e-03 3.79537171e-04 4.34335669e-04 1.46176111e-03
 1.81276561e-03 6.18689578e-04 3.64098379e-04 2.99725974e-04
 4.41121691e-02 8.69079731e-04 4.63531291e-04 4.38290542e-04
 3.95782673e-04 2.94414239e-04 1.45253916e-03 4.50655418e-04
 2.88441108e-04 3.87161378e-04 3.66006522e-04 9.21210440e-04
 9.62048419e-04 3.31125955e-03 4.62450237e-04 6.10613510e-04
 2.83632763e-04 3.13751494e-04 5.73204750e-04 9.66607950e-04
 3.66662549e-04 1.86712979e-03 2.70711139e-04 3.30637939e-04
 4.16621222e-04 4.80582597e-03 2.74162517e-04 8.27324123e-04
 4.37910085e-03 2.99351426e-04 3.04995191e-03 4.56892388e-04
 4.52613038e-04 5.36569442e-04 5.52969777e-04 6.79600369e-04
 3.71211781e-04 2.55522024e-04 4.41834960e-04 4.99009844e-04
 3.31654890e-04 2.71205205e-03 3.09116327e-04 3.83554603e-04
 6.33114795e-04 8.27489363e-04 3.20148895e-03 4.79803828e-03
 4.39460248e-04 4.14850296e-04 2.98175198e-04 2.74712755e-04
 3.68523367e-04 2.96752158e-04 3.48303317e-04 4.35518144e-04
 9.34146132e-04 3.58313083e-04 3.62306067e-03 5.94764207e-04
 6.82957434e-03 5.11217980e-04 3.57485411e-03 1.82299160e-03
 1.36522236e-03 2.98418005e-04 3.56942147e-04 4.66868862e-04
 2.08039654e-03 4.00715660e-04 3.05312182e-04 6.50790027e-04
 3.29660141e-04 5.15219091e-04 4.36037488e-03 9.54040651e-04
 1.30418181e-03 6.11263987e-04 1.46496999e-03 6.57213027e-03
 6.61310109e-04 1.45812009e-03 3.38598482e-04 4.34048838e-04
 9.71587663e-04 6.17032432e-04 6.01449533e-04 2.82616984e-04
 1.31176835e-03 2.20882394e-03 2.90881325e-04 3.70669677e-04
 5.27466076e-03 3.82776893e-04 3.10523565e-04 2.08058578e-03
 2.09262135e-03 2.48706180e-04 5.54920785e-04 5.74331637e-04
 2.73627368e-04 6.16772075e-04 4.97419818e-04 6.56238933e-04
 1.77115720e-03 4.01074541e-04 3.94210327e-04 5.49362358e-04
 3.58439557e-04 3.61686777e-04 2.54864542e-04 3.61748320e-04
 3.27643619e-04 1.04546293e-03 2.83205488e-04 3.29696244e-04
 2.67081383e-04 7.56376112e-04 2.45446931e-04 3.36356162e-04
 5.80408253e-04 4.76466196e-04 4.74342320e-04 2.65400565e-04
 8.40893288e-03 6.85550045e-04 4.49351575e-04 1.48871279e-03
 4.62762831e-04 4.13256732e-04 8.10369204e-04 3.23475823e-04
 4.19468061e-04 2.99378497e-04 7.42618687e-04 2.47048695e-04
 3.57538872e-04 2.65767072e-04 3.98523666e-04 4.16083661e-04
 1.19607817e-03 3.06136805e-04 7.32927370e-04 1.50379308e-03
 2.62762594e-04 4.53491898e-04 3.78522623e-04 4.31704608e-04
 1.75004759e-03 1.25032776e-03 3.27301325e-04 3.12266393e-04
 7.51865473e-04 4.88383039e-04 3.13137587e-04 3.88165212e-04
 3.38210593e-04 1.70122071e-03 2.52645185e-04 3.84139581e-04
 3.15484246e-04 4.05445951e-03 5.90037902e-04 5.01378207e-03
 2.75757842e-04 3.74035113e-04 1.35626379e-03 4.46171104e-03
 5.61642764e-04 4.69177748e-04 3.63259146e-04 3.00602698e-03
 4.98867073e-04 8.40009249e-04 5.37346944e-04 4.74872303e-03
 3.07012536e-04 4.71755072e-03 2.79907198e-04 3.72254610e-04
 5.46859136e-04 7.01695188e-04 3.62247955e-04 2.62301041e-04
 2.40390396e-03 3.82269562e-04 4.59987440e-03 4.79291848e-04
 3.23994832e-04 8.71842545e-04 7.31086891e-04 2.91781864e-04
 2.53599743e-04 3.79510317e-04 5.19125563e-04 2.00023531e-03
 3.88562380e-04 3.87007973e-03 6.70528630e-04 7.00383138e-04
 3.11933556e-04 3.30141785e-04 7.90048123e-04 1.51745140e-03
 6.73859080e-04 5.51368068e-04 3.40206222e-01 4.09083815e-03
 3.76891564e-04 4.40940057e-04 6.65542584e-04 5.19568296e-04
 1.39848688e-03 4.80237782e-04 2.45207300e-04 5.29121233e-03
 2.64659410e-04 1.55945230e-03 5.47967301e-04 7.98754034e-04
 8.95172430e-04 2.07589838e-03 3.09065176e-04 3.57159833e-04
 4.11161736e-04 1.14456899e-03 2.32668001e-03 3.03037007e-04
 4.16294004e-04 2.17498701e-03 3.06452930e-04 2.15205450e-02
 2.58303543e-04 6.51516359e-04 6.42360248e-04 5.02449866e-04
 8.80919004e-04 3.89958234e-04 1.56892155e-03 3.20696952e-04
 8.40111749e-03 3.62088526e-04 8.85107059e-03 3.09165330e-04
 4.37383808e-04 7.56985962e-04 5.16424609e-03 7.23450017e-03
 3.30088735e-04 2.72886104e-04 1.36296752e-03 2.53349179e-04
 7.25696693e-04 3.72832272e-03 4.24933142e-04 4.78915688e-04
 2.88675846e-04 3.51261179e-04 1.90341225e-03 4.30589697e-04
 2.55641258e-04 3.28793025e-04 7.57311788e-04 2.28638966e-03
 3.68079850e-04 4.08505951e-04 5.71515284e-04 4.47527006e-04
 4.77598892e-04 5.20534070e-04 2.48821737e-04 6.67810160e-04
 2.80479901e-04 4.77746527e-04 2.58099360e-04 5.27701037e-04
 1.01523059e-03 9.41717491e-03 1.76423553e-03 4.64299667e-04
 4.18435904e-04 4.23968342e-04 4.18389898e-03 2.93162192e-04
 1.30443196e-03 6.57153863e-04 4.94020048e-04 2.73543728e-04
 9.55453033e-03 3.55310572e-04 1.66307074e-03 3.00990435e-03
 3.59473297e-04 1.14708882e-03 3.37876852e-04 2.43977525e-03
 7.33451538e-04 4.82836805e-04 7.05398513e-04 1.03345980e-03
 1.28729542e-03 4.88666974e-03 8.72852756e-04 2.70601299e-04
 7.63579349e-04 7.04879953e-04 1.70529186e-03 2.98582889e-04
 8.80295037e-04 9.45019235e-04 2.45798045e-04 6.77037524e-04
 3.80978771e-04 3.70850286e-04 3.85126597e-04 2.60652287e-04
 3.03591296e-04 3.83213722e-04 2.42449874e-03 2.62531118e-04
 8.73664514e-04 5.34530670e-04 6.26206579e-04 2.10418500e-03
 1.13529876e-03 3.99757686e-04 4.80677677e-03 2.64539355e-04
 3.57468005e-04 8.71438026e-04 3.05030313e-04 3.38028708e-04
 3.24752940e-04 2.88237398e-04 6.36844746e-04 6.11920957e-04
 3.09884554e-04 3.78816348e-03 3.18898689e-03 3.90350035e-04
 1.21977609e-03 4.34112472e-04 7.11096691e-04 3.65261080e-04
 3.18131957e-03 5.19152658e-04 2.25585445e-03 1.31315613e-03
 6.28507762e-04 2.60907670e-04 2.81780237e-04 4.40587385e-04
 9.29770417e-04 6.18218125e-04 5.46538127e-04 4.73675710e-04
 3.33895448e-03 2.96708169e-04 5.68386272e-04 2.79242822e-04
 9.32904360e-04 9.40540213e-04 8.73715661e-04 1.08553899e-03
 3.69807333e-04 4.68858814e-04 2.57848006e-04 1.04619269e-03
 2.96015914e-04 5.94117120e-04 2.48185127e-04 2.85245992e-04
 3.28402505e-04 1.18616675e-03 9.29303308e+00 5.49886390e-04
 4.53443430e-04 8.98135826e-04 8.49355748e-04 3.41861254e-04
 8.02769323e-04 2.59610601e-04 6.57754109e-04 3.04390099e-04
 3.65636458e-04 1.00797425e-03 2.89614985e-04 7.52516897e-04
 3.52428481e-04 2.57021878e-04 3.45776714e-04 4.25440996e-04
 4.95424860e-04 4.93146020e-04 6.31607100e-04 3.37135500e-04
 6.85598420e-04 3.41713065e-04 3.47394867e-03 4.75862673e-04
 6.99887261e-04 5.28215197e-04 1.94112686e+00 2.43637158e-03
 3.97398536e-04 6.39634071e-04 2.58069091e-04 3.64716610e-04
 3.04666472e-04 3.11696796e-04 4.26392216e-04 2.86567595e-04
 4.63680529e-04 3.52411077e-04 3.13099501e-04 2.79049787e-04
 3.02202980e-04 2.49901590e-04 3.82921309e-03 3.37371691e-04
 8.91423632e-04 5.39941052e-03 4.48885326e-04 4.07443246e-04
 7.33543785e-04 4.67380034e-04 5.37602071e-03 4.11588986e-04
 6.19731351e-04 4.26044123e-03 5.08720585e-04 1.05689846e-03
 3.72293013e-04 4.70739582e-04 6.20712803e-04 2.54363847e-04
 3.82760608e-04 1.71267508e-03 6.26260569e-04 5.61326092e-04
 3.15994858e-04 7.27738337e-04 2.63183174e-04 4.69703159e-04
 4.54973559e-04 2.47869166e-04 2.98776045e-04 3.37824715e-04
 6.52307238e-04 2.99002790e-04 8.14543684e-04 1.14853992e-03
 3.06111547e-04 2.70671448e-04 3.58912851e-03 2.95878872e-03
 2.70420200e-04 3.63572997e-04 5.58891976e-04 3.11332134e-04
 4.47067682e-04 2.66852940e-04 2.49295725e-04 4.56368582e-04
 4.03073713e-04 7.79036651e-04 4.45589887e-04 2.63499974e-04
 2.71455846e-04 3.58633169e-04 6.22896251e-04 3.41261030e-03
 5.59483680e-04 2.79288693e-04 3.28092543e-04 4.93296139e-04
 6.18207619e-04 3.58118793e-04 9.22919457e-03 1.01698828e-03
 4.46036129e-03 2.69571827e-04 2.69413324e-04 6.09778570e-04
 4.16621065e-04 3.27362048e-04 7.84633461e-04 5.69964234e-04
 2.94273342e-04 7.88776081e-04 5.73427166e-04 7.69329874e-04
 5.63132132e-04 3.82785182e-04 4.80224041e-04 4.67819955e-03
 2.93183311e-04 2.85653497e-04 5.26704673e-04 4.49319583e-03
 2.60544020e-04 3.03403553e-04 3.83932886e-04 2.73482078e-04
 6.36864469e-04 6.29663818e-03 1.93829988e-03 2.65873474e-04
 4.29386083e-04 7.98255073e-04 3.45317099e-04 2.92108251e-04
 6.85901435e-04 2.00085431e-03 2.63800696e-04 3.02523732e-04
 3.95086387e-04 9.88739159e-04 9.79618143e-03 2.65917969e-04
 2.95889252e-04 3.34933476e-04 3.33607404e-04 4.39227563e-03
 3.37565380e-04 2.84080493e-04 2.57103220e-04 4.87948003e-04
 2.74950736e-04 7.35561526e-04 3.28153321e-04 3.55074314e-04
 3.34911661e-04 4.86693557e-04 3.94888958e-04 4.63609727e-04
 3.39861959e-04 3.25750411e-04 3.37174171e-04 3.59227491e-04
 4.73408109e-04 1.00420394e-03 3.34189302e-04 4.73270827e-04
 5.01370224e-04 3.21689031e-04 4.10616719e-03 2.63910823e-03
 3.82457950e-04 1.13528053e-03 1.02591202e-03 2.50920965e-04
 2.49476617e-04 3.74311851e-04 5.73641592e-04 4.62242764e-04
 3.82551814e-03 6.97545209e-04 4.25620530e-03 4.76380839e-04
 4.60194360e-04 2.68900001e-04 4.65780634e-04 7.61949117e-04
 5.74783284e-04 2.93135921e-04 3.79753843e-04 7.05919878e-04
 3.92907105e-04 2.64060448e-04 5.82742792e-04 3.13448748e-04
 3.92214076e-04 6.46059214e-04 4.25092730e-04 4.81368983e-04
 7.61617768e-04 3.74740040e-04 2.69641578e-04 3.59004818e-04
 4.92355790e-04 3.71278909e-04 6.40900628e-03 3.12586210e-04
 2.58056383e-04 3.11797869e-03 4.47675955e-04 7.82446476e-04
 2.59928280e-04 2.50198218e-04 5.57714432e-04 3.13549223e-04
 3.45425403e-04 3.04159822e-04 6.08585190e-04 4.45880833e-04
 5.86878759e-03 6.92692489e-04 9.23881597e-01 5.13395776e-04
 3.02779635e-03 3.41076341e-04 5.53810879e-04 9.13558594e-04
 3.36329717e-04 1.48721106e-03 1.06533957e-03 4.87951157e-04
 4.78071534e-04 2.82509696e-04 2.84646006e-03 2.34202207e-03
 4.97507204e-04 3.66535812e-04 2.70075220e-04 1.09193708e-03
 2.67826847e-04 3.65926883e-04 2.71532543e-04 3.79689580e-04
 2.95138559e-04 6.74076825e-04 3.81687946e-04 3.67945522e-04
 1.43888506e-03 3.22455718e-04 3.78351328e-04 3.17660904e-04
 4.11798201e-04 3.06942337e-04 9.25289483e-04 2.76054872e-04
 3.44925121e-04 9.70486689e-04 5.76658032e-03]

Auch in Modell lm2 werden mögliche outliers und Punkte mit high leverage identifiziert und aus dem Datenset entfernt.

#drop outliers identified by Cook's distance
train_dataset2 = train_dataset.drop(train_dataset.index[out_d])
train_dataset2.info()
<class 'pandas.core.frame.DataFrame'>
Int64Index: 15484 entries, 14185 to 17879
Data columns (total 12 columns):
 #   Column                Non-Null Count  Dtype   
---  ------                --------------  -----   
 0   housing_median_age    15484 non-null  int64   
 1   total_rooms           15484 non-null  int64   
 2   total_bedrooms        15484 non-null  int64   
 3   population            15484 non-null  int64   
 4   households            15484 non-null  int64   
 5   median_income         15484 non-null  float64 
 6   median_house_value    15484 non-null  float64 
 7   ocean_proximity       15484 non-null  category
 8   people_per_household  15484 non-null  float64 
 9   rooms_per_household   15484 non-null  float64 
 10  bedrooms_per_people   15484 non-null  float64 
 11  bedrooms_per_rooms    15484 non-null  float64 
dtypes: category(1), float64(6), int64(5)
memory usage: 1.4 MB
# Fit Model
lm2 = smf.ols(formula='median_house_value ~ median_income + ocean_proximity + bedrooms_per_rooms + people_per_household + housing_median_age', data=train_dataset2).fit()
lm2.summary()
OLS Regression Results
Dep. Variable: median_house_value R-squared: 0.746
Model: OLS Adj. R-squared: 0.746
Method: Least Squares F-statistic: 9085.
Date: Wed, 19 Jan 2022 Prob (F-statistic): 0.00
Time: 14:04:48 Log-Likelihood: -1.9034e+05
No. Observations: 15484 AIC: 3.807e+05
Df Residuals: 15478 BIC: 3.807e+05
Df Model: 5
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept -4.951e+04 4086.781 -12.115 0.000 -5.75e+04 -4.15e+04
ocean_proximity[T.NEAR COAST] 5.763e+04 1060.730 54.330 0.000 5.56e+04 5.97e+04
median_income 4.857e+04 368.630 131.763 0.000 4.78e+04 4.93e+04
bedrooms_per_rooms 3.637e+05 1.19e+04 30.623 0.000 3.4e+05 3.87e+05
people_per_household -2.752e+04 578.566 -47.572 0.000 -2.87e+04 -2.64e+04
housing_median_age 923.2678 36.548 25.262 0.000 851.629 994.907
Omnibus: 1031.972 Durbin-Watson: 1.988
Prob(Omnibus): 0.000 Jarque-Bera (JB): 1351.257
Skew: 0.610 Prob(JB): 3.79e-294
Kurtosis: 3.778 Cond. No. 912.


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

Die Modellperformance liegt nun bei 74,6% und ist somit bisher der beste erzielte Wert der Kennzahl R-squared.
F-statistics hat sich ebenfalls verbessert.
Alle anderen Werte liegen in einem ähnlichen Wertebereich.

fig = sm.graphics.plot_partregress_grid(lm2)
fig.tight_layout(pad=0.2)
_images/regression_128_0.png

Alle Features zeigen einen Einfluss auf median_house_value.

fig = sm.graphics.plot_regress_exog(lm2, "housing_median_age")
fig.tight_layout(pad=1.0)
_images/regression_130_0.png

Das Partial Regression Plot und CCPR Plot zeigen eine geringe lineare Korrelation von housing_median_age in Verbindung mit den anderen unabhängigen Variablen auf median_house_value.
Das Residual Plot zeigt eine deutliche Homoskedastizität.

fig = sm.graphics.plot_regress_exog(lm2, "ocean_proximity[T.NEAR COAST]")
fig.tight_layout(pad=1.0)
_images/regression_132_0.png

Durch die Zusammenlegung der Ausprägungen <1H OCEAN und NEAR COST hat ds Feature nun eine positive lineare Abhängigkeit zu median_house_value. Die Fehlerstreuung ist bei Distrikten in Meerehöhe etwas größer. Auch hier wird vorsichtshalber der Breusch-Pagan Lagrange Multiplier Test zum Ausschluss von Heteroskedastizität durchgeführt.

name = ['Lagrange multiplier statistic', 'p-value', 'f-value', 'f p-value']
test = sm.stats.het_breuschpagan(lm2.resid, lm2.model.exog)
lzip(name, test)
[('Lagrange multiplier statistic', 1095.2984925363032),
 ('p-value', 1.3939938205969558e-234),
 ('f-value', 235.64364107050312),
 ('f p-value', 2.5613416182330695e-243)]

Da alle p-values über 0,05 liegen, kann Heteroskedastizität ausgeschlossen werden.

In der linearen Regression muss eine Colinearität der unabhängigen Variablen unbedingt vermieden werden. Mithilfe des Variance Inflation Factors (VIF) kann der Grad der Colinearität der im Modell verwendeten Variablen berechnet werden. Der kleinstmögliche Wert ist 1. Problematisch wird es ab Werten > 5.

# choose features and add constant
features = add_constant(train_dataset2[['median_income', 'bedrooms_per_rooms', 'people_per_household', 'housing_median_age']])
# create empty DataFrame
vif = pd.DataFrame()
# calculate vif
vif["VIF Factor"] = [variance_inflation_factor(features.values, i) for i in range(features.shape[1])]
# add feature names
vif["Feature"] = features.columns

vif.round(2)
C:\ProgramData\Anaconda3\lib\site-packages\statsmodels\tsa\tsatools.py:142: FutureWarning:

In a future version of pandas all arguments of concat except for the argument 'objs' will be keyword-only
VIF Factor Feature
0 83.61 const
1 1.80 median_income
2 1.79 bedrooms_per_rooms
3 1.01 people_per_household
4 1.03 housing_median_age

Der Wert der Konstanten kann ignoriert werden. Ansonsten ist die Colinearität der Variablen eher gering.

Fazit zu Modellen 1 und 2:
Die Modelle lm und lm2 sind sehr ähnlich. R-squared ist bei lm2 ca. 1% besser, auch der Wert von F-statistics hat sich von lm zu lm2 verbessert. AIC und BIC haben sich dagegen minimal verschlechtert Das Partial Regression Plot der Variablen median_housing_age zeigt zwar nur eine geringe Linearität, dennoch scheint die Aufnahme der Variablen das Modell etwas zu verbessern. Daher verwenden die nachfolgenden Modelle die in lm2 verwendeten Features.

Model 2 - Lasso Regression

Nachfolgend wird eine Lasso Regression mit statsmodels fit_regularized() Methode trainiert. Der Wert von alpha wird auf 5 gesetzt.

lasso1 = smf.ols(formula='median_house_value ~ median_income + ocean_proximity + bedrooms_per_rooms + people_per_household + housing_median_age', data=train_dataset2).fit_regularized(method='sqrt_lasso', alpha=5.0, L1_wt=1.0, max_iter=1000)
print("Parameters: ", lasso1.params)
Parameters:  Intercept                        -44272.058534
ocean_proximity[T.NEAR COAST]     58127.019057
median_income                     48170.671059
bedrooms_per_rooms               346772.197084
people_per_household             -27621.154742
housing_median_age                  917.516632
dtype: float64

Da die Summary() Methode für statsmodels fit_regularized() Medelle nicht verfügbar ist, müssen wichtige Kennzahlen selbst berechnet werden.

# calculate the mean 
train_dataset2['average'] = train_dataset2["median_house_value"].mean()

# calculate error (observation - average) and assign it to dataframe
train_dataset2 = train_dataset2.assign(error = (train_dataset2['median_house_value'] - train_dataset2['average']))

# calculate squared error and assign it to dataframe
train_dataset2 = train_dataset2.assign(error_sq = (train_dataset2['median_house_value'] - train_dataset2['average'])**2)


# obtain predicted value from model lasso
train_dataset2['pred'] = lasso1.predict()

# obtain the residuals from statsmodel (resid)
train_dataset2['error_2'] = train_dataset2['median_house_value'] - train_dataset2['pred']

# square the residuals 
train_dataset2['error_sq_2'] = train_dataset2['error_2']**2

# show df
train_dataset2.head(5)
housing_median_age total_rooms total_bedrooms population households median_income median_house_value ocean_proximity people_per_household rooms_per_household bedrooms_per_people bedrooms_per_rooms average error error_sq pred error_2 error_sq_2
14185 37 2176 418 1301 375 2.8750 98900.0 NEAR COAST 3.469333 5.802667 0.321291 0.192096 196573.25704 -97673.25704 9.540065e+09 157080.171533 -58180.171533 3.384932e+09
6125 20 3158 684 2396 713 3.5250 153000.0 NEAR COAST 3.360449 4.429173 0.285476 0.216593 196573.25704 -43573.25704 1.898629e+09 184295.786398 -31295.786398 9.794262e+08
14095 11 2393 726 1905 711 1.3448 91300.0 NEAR COAST 2.679325 3.365682 0.381102 0.303385 196573.25704 -105273.25704 1.108246e+10 119926.953213 -28626.953213 8.195025e+08
14359 52 1260 202 555 209 7.2758 345200.0 NEAR COAST 2.655502 6.028708 0.363964 0.160317 196573.25704 148626.74296 2.208991e+10 394291.589328 -49091.589328 2.409984e+09
18004 32 2930 481 1336 481 6.4631 344100.0 NEAR COAST 2.777547 6.091476 0.360030 0.164164 196573.25704 147526.74296 2.176414e+10 334755.756942 9344.243058 8.731488e+07
# calculate sum of squared error (which is in case of the mean the total error)
TSS = train_dataset2.error_sq.sum()
# print output
print('Sum of squared error (TSS) of model 1:', TSS)

# Sum of squared residuals (SS_R)
SSR = train_dataset2['error_sq_2'].sum()
print('Sum of squared residuals (SSR) of model 1:',SSR)
Sum of squared error (TSS) of model 1: 169695617828075.0
Sum of squared residuals (SSR) of model 1: 43134248561958.78
# Explained sum of squares  (SS_M = TSS - SS_R)
SSM = TSS - SSR
print('SSM of Lasso Model: ', SSM)

# R_Squared: explained sum of squared residuals
R_squared = SSM / TSS
print('R squared Lasso Model:', R_squared)
SSM of Lasso Model:  126561369266116.22
R squared Lasso Model: 0.7458140103201739

R-squared ist durch die Regularisierung nochmal etwas besser geworden.
Für einen optimalen Wert von alpha müsste eigentlich eine Cross-Validation durchgeführt werden. Darauf wird an dieser Stelle jedoch verzichtet.

Model 3 - Cubic Spline

Nun wird eine Cubic Spline unter Verwendung von Patsy und Statsmodels trainiert. Eine Verwendung von mehreren unabhängigen Variablen ist dabei nicht ohne Weiteres möglich. Daher wird nachfolgend lediglich die Variable mit der höchsten Korrelation zu median_house_value verwendet. Dabei handelt es sich um median_income.

X_train = train_dataset2['median_income']
y_train = train_dataset2['median_house_value']
X_test = test_dataset['median_income']
y_test = test_dataset['median_house_value']
# Generating cubic spline with 3 knots at 25, 40 and 60
transformed_x = dmatrix(
            "bs(train, knots=(3,6,9), degree=3, include_intercept=False)", 
                {"train": X_train},return_type='dataframe')
# Fitting generalised linear model on transformed dataset
spline1 = sm.GLM(y_train, transformed_x).fit()
spline1.predict()
spline1.summary()
Generalized Linear Model Regression Results
Dep. Variable: median_house_value No. Observations: 15484
Model: GLM Df Residuals: 15477
Model Family: Gaussian Df Model: 6
Link Function: identity Scale: 4.5931e+09
Method: IRLS Log-Likelihood: -1.9421e+05
Date: Wed, 19 Jan 2022 Deviance: 7.1087e+13
Time: 14:14:30 Pearson chi2: 7.11e+13
No. Iterations: 3
Covariance Type: nonrobust
coef std err z P>|z| [0.025 0.975]
Intercept 1.08e+05 8224.837 13.135 0.000 9.19e+04 1.24e+05
bs(train, knots=(3, 6, 9), degree=3, include_intercept=False)[0] -3.733e+04 1.09e+04 -3.439 0.001 -5.86e+04 -1.61e+04
bs(train, knots=(3, 6, 9), degree=3, include_intercept=False)[1] 7.416e+04 7387.742 10.038 0.000 5.97e+04 8.86e+04
bs(train, knots=(3, 6, 9), degree=3, include_intercept=False)[2] 1.573e+05 9964.123 15.782 0.000 1.38e+05 1.77e+05
bs(train, knots=(3, 6, 9), degree=3, include_intercept=False)[3] 4.465e+05 1.06e+04 42.023 0.000 4.26e+05 4.67e+05
bs(train, knots=(3, 6, 9), degree=3, include_intercept=False)[4] 3.512e+05 1.97e+04 17.846 0.000 3.13e+05 3.9e+05
bs(train, knots=(3, 6, 9), degree=3, include_intercept=False)[5] 4.248e+05 4.12e+04 10.307 0.000 3.44e+05 5.06e+05
# obtain predicted value from model spline
train_dataset2['pred_spline'] = spline1.predict()

# obtain the residuals from statsmodel (resid)
train_dataset2['error_2_spline'] = train_dataset2['median_house_value'] - train_dataset2['pred_spline']

# square the residuals 
train_dataset2['error_sq_2_spline'] = train_dataset2['error_2_spline']**2
housing_median_age total_rooms total_bedrooms population households median_income median_house_value ocean_proximity people_per_household rooms_per_household ... bedrooms_per_rooms average error error_sq pred error_2 error_sq_2 pred_spline error_2_spline error_sq_2_spline
14185 37 2176 418 1301 375 2.8750 98900.0 NEAR COAST 3.469333 5.802667 ... 0.192096 196573.25704 -97673.25704 9.540065e+09 157080.171533 -58180.171533 3.384932e+09 154242.775962 -55342.775962 3.062823e+09
6125 20 3158 684 2396 713 3.5250 153000.0 NEAR COAST 3.360449 4.429173 ... 0.216593 196573.25704 -43573.25704 1.898629e+09 184295.786398 -31295.786398 9.794262e+08 183202.797434 -30202.797434 9.122090e+08
14095 11 2393 726 1905 711 1.3448 91300.0 NEAR COAST 2.679325 3.365682 ... 0.303385 196573.25704 -105273.25704 1.108246e+10 119926.953213 -28626.953213 8.195025e+08 96483.272131 -5183.272131 2.686631e+07
14359 52 1260 202 555 209 7.2758 345200.0 NEAR COAST 2.655502 6.028708 ... 0.160317 196573.25704 148626.74296 2.208991e+10 394291.589328 -49091.589328 2.409984e+09 372203.127060 -27003.127060 7.291689e+08
18004 32 2930 481 1336 481 6.4631 344100.0 NEAR COAST 2.777547 6.091476 ... 0.164164 196573.25704 147526.74296 2.176414e+10 334755.756942 9344.243058 8.731488e+07 319584.526302 24515.473698 6.010085e+08

5 rows × 21 columns

# Sum of squared residuals (SSR)
SSR_spline = train_dataset2['error_sq_2_spline'].sum()
print('Sum of squared residuals (SSR) of Cubic Spline:',SSR_spline)

# Explained sum of squares  (SS_M = TSS - SS_R)
SSM_spline = TSS - SSR_spline
print('SSM of Cubic Spline: ', SSM_spline)

# R_Squared: explained sum of squared residuals
R_squared_spline = SSM_spline / TSS
print('R squared of Cubic Spline:', R_squared_spline)
Sum of squared residuals (SSR) of Cubic Spline: 71087433258433.94
SSM of Cubic Spline:  98608184569641.06
R squared of Cubic Spline: 0.581088573952126

R squared ist schlechter geworden als in den vorherigen Modellen. Allerdings liegt dies darin begründet, dass nur ein Feature verwendet wurde. Ein direkter Vergleich ist also nicht möglich.

# Create observations
xp = np.linspace(X_test.min(),X_test.max(), 100)
# Make some predictions
pred = spline1.predict(dmatrix("bs(xp, knots=(3, 6, 9), include_intercept=False)", {"xp": xp}, return_type='dataframe'))

# plot
sns.scatterplot(x=X_train, y=y_train)

plt.plot(xp, pred, label='Cubic spline with degree=3 (3 knots)', color='orange')
plt.legend();
_images/regression_158_0.png

Model 4 - Natural Spline

transformed_x3 = dmatrix("cr(train,df = 3)", {"train": X_train}, return_type='dataframe')

spline2 = sm.GLM(y_train, transformed_x3).fit()
# obtain predicted value from model spline
train_dataset2['pred_spline'] = spline2.predict()

# obtain the residuals from statsmodel (resid)
train_dataset2['error_2_spline'] = train_dataset2['median_house_value'] - train_dataset2['pred_spline']

# square the residuals 
train_dataset2['error_sq_2_spline'] = train_dataset2['error_2_spline']**2
# Sum of squared residuals (SSR)
SSR_spline = train_dataset2['error_sq_2_spline'].sum()
print('Sum of squared residuals (SSR) of Natural Spline:',SSR_spline)

# Explained sum of squares  (SS_M = TSS - SS_R)
SSM_spline = TSS - SSR_spline
print('SSM of Natural Spline: ', SSM_spline)

# R_Squared: explained sum of squared residuals
R_squared_spline = SSM_spline / TSS
print('R squared of Natural Spline:', R_squared_spline)
Sum of squared residuals (SSR) of Natural Spline: 72042025549993.14
SSM of Natural Spline:  97653592278081.86
R squared of Natural Spline: 0.5754632531349064

R-squared ist im Vergleich zur Cubic Spline nochmal gesunken. Auch hier wird im Modell nur eine unabhängige Variable verwendet.

# Make predictions
pred = spline2.predict(dmatrix("cr(xp, df=3)", {"xp": xp}, return_type='dataframe'))

# plot
sns.scatterplot(x=X_train, y=y_train)
plt.plot(xp, pred, color='orange', label='Natural spline with df=3')
plt.legend()
<matplotlib.legend.Legend at 0x222f091b2b0>
_images/regression_164_1.png

Scikit Learn

Data Preprocessing Pipeline

Initialisierung der Preprocessing Pipeline. Mithilfe einer Preprocessing Pipelines müssen die Daten nicht manuell transformiert werden, da diese Aufgaben von der Pipeline übernommen werden.
Für numerische Variablen ist dies beispielsweise die Normalisierung der Variablen, sowie die Belegung nicht befüllter Datensätze mit dem Durchschnittswert des Features.
Bei kategorialen Variablen werden fehlnde Werte durch eine Konstante belegt und automatisiert Dummy Variablen für die verschiedenen Ausprägungen der Variablen erzeugt.

# for numeric features
numeric_transformer = Pipeline(steps=[
    ('imputer', SimpleImputer(strategy='median')),
    ('scaler', StandardScaler()) #Normalisierung (X=(X-Mittelwert) / Standardabweichung)
    ])
# for categorical features  
categorical_transformer = Pipeline(steps=[
    ('imputer', SimpleImputer(strategy='constant', fill_value='missing')), #Konstante bei fehlenden Werten reinschreiben
    ('onehot', OneHotEncoder(handle_unknown='ignore')) #pro Featureausprägung eine Spalte, zutreffendes Feature hat dann Wert 1 und die anderen Wert 0
    ])
# Pipeline
preprocessor = ColumnTransformer(transformers=[
    ('num', numeric_transformer, selector(dtype_exclude="category")),
    ('cat', categorical_transformer, selector(dtype_include="category"))
        ])

Datenvorbereitung

Für Scikit Learn muss die abhängige Variable von den unabhängigen getrennt werden. Als Datenbasis wird hier auf train_dataset2 zugegriffen. In diesem Datenset wurden die Outlier nach Cook’s Distance bereits entfernt. Als Features sollen die durch die Filter Methode identifizierten und bereits in statsmodels verwendeten Features genutzt werden.

# create label y
y = train_dataset2['median_house_value']

# create features
X= train_dataset2[['median_income', 'ocean_proximity', 'bedrooms_per_rooms', 'people_per_household', 'housing_median_age']]

Scikit Learn bietet eine eigene Methode zum Datensplitting. Das Datenset train_dataset2 aus statsmodels wird nun nochmal in Test- und Trainingsdaten aufgeteilt. Dadurch ist eine Validierung der Modelle bereits mit diesem Testset möglich, bevor am Ende dieses Notebooks mit den eigentlichen Testdaten das beste Modell evaluiert wird.

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=10)
print("Training size:", len(X_train))
print("Test size:", len(X_test))
Training size: 12387
Test size: 3097

Model 5 - Simple Regression

Nun wird ein Pipeline Object erzeugt, welches zunächst die Data Preprocessing Pipeline ausführt und anschließend das Regressionsmodell.

# Create pipeline with model
lm_pipe = Pipeline(steps=[
    ('preprocessor', preprocessor),
    ('lm', LinearRegression())
                        ])

# Fit model
lm_pipe.fit(X_train, y_train)
Pipeline(steps=[('preprocessor',
                 ColumnTransformer(transformers=[('num',
                                                  Pipeline(steps=[('imputer',
                                                                   SimpleImputer(strategy='median')),
                                                                  ('scaler',
                                                                   StandardScaler())]),
                                                  <sklearn.compose._column_transformer.make_column_selector object at 0x00000222F09495B0>),
                                                 ('cat',
                                                  Pipeline(steps=[('imputer',
                                                                   SimpleImputer(fill_value='missing',
                                                                                 strategy='constant')),
                                                                  ('onehot',
                                                                   OneHotEncoder(handle_unknown='ignore'))]),
                                                  <sklearn.compose._column_transformer.make_column_selector object at 0x00000222F0949940>)])),
                ('lm', LinearRegression())])

In der Variable features_names werden die Namen der Features gespeichert.

feature_names = np.concatenate((X_train.drop(['ocean_proximity'], axis=1).columns.to_numpy(), lm_pipe.named_steps['preprocessor'].transformers_[1][1]['onehot'].get_feature_names_out()))
feature_names
array(['median_income', 'bedrooms_per_rooms', 'people_per_household',
       'housing_median_age', 'x0_INLAND', 'x0_NEAR COAST'], dtype=object)

Der eben erzeugte feature_names Array kann nun mit den Koeffizienten des Modells ergänzt werden, um die Koeffizienten je Feature auszugeben.

# Obtain model coefficients and names
print('Koeffizienten und Namen: \n', list(zip(lm_pipe.named_steps['lm'].coef_, feature_names)))
Koeffizienten und Namen: 
 [(84852.55471482786, 'median_income'), (18735.049763174888, 'bedrooms_per_rooms'), (-19963.944594294095, 'people_per_household'), (11488.308364831402, 'housing_median_age'), (-28587.94998804061, 'x0_INLAND'), (28587.949988040604, 'x0_NEAR COAST')]

Das trainierte Modell wird nun für die Berechnung der median_house_values der Trainingsdaten verwendet.

y_pred = lm_pipe.predict(X_train)
y_pred
array([113004.32375212, 163416.09509259, 104877.65553309, ...,
        45560.80974401, 306094.8367317 , 117268.25027873])
# Metrics trainings set
print('R squared training set:', round(skm.r2_score(y_train, y_pred)*100, 2))
print('MSE training set:', round(skm.mean_squared_error(y_train, y_pred), 2))
print('RMSE training set:', round(skm.mean_squared_error(y_train, y_pred, squared=False), 2))
R squared training set: 74.47
MSE training set: 2777618420.95
RMSE training set: 52703.12

Das Modell hat einen ähnlichen R-squared wie die OLS Regression in statsmodels. Der Root Mean Squared Error des Modells beträgt 52703,12$. Das bedeutet, dass der eigentliche Wert von median_house_value im Schnitt um diesen Betrag von dem durch das Modell berechneten Ergebnis abweicht.

Im Folgenden wird das Modell mit den im Scikit Learn Datensplit erzeugten Validierungsdaten getestet.

# Validation data
pred = lm_pipe.predict(X_test)
r2_test = round(skm.r2_score(y_test, pred)*100, 2)
print('R squared validation set:', r2_test)
mse_test =skm.mean_squared_error(y_test, pred)
print('MSE validation set', round(mse_test, 2))
print('RMSE training set:', round(skm.mean_squared_error(y_test, pred, squared=False), 2))
R squared validation set: 75.04
MSE validation set 2816769494.14
RMSE training set: 53073.25

Der R-squared Wert hat sich mit den Validierungsdaten verbessert. Dies ist sehr ungewöhnliche, da sich unter Verwendung unbekannter Variablen die Kennzahlen normalerweise etwas verschlechtern. Allerdings wurde mit diesen Validierungsdaten bereits EDA, Feature Engineering und auch Cook’s Distance Outlier Entfernung durchgeführt. Die Daten sind somit nicht wirklich unbekannt.
Der durchschnittliche Fehler MSE und auch die durchschnittliche Abweichung von berechneten zu tatsächlichen Modellen sind etwas schlechter geworden. Allerdings nur geringfügig.

Model 6 - Lasso

vorgegebenes Alpha

In diesem Unterkapitel wird eine Lasso Regression für das Modell durchgeführt. Der Wert für alpha wird initial mit 1 belegt.

# Create pipeline with model
lasso_pipe = Pipeline(steps=[
    ('preprocessor', preprocessor),
    ('lasso', Lasso(alpha=1, max_iter=10000))
                        ])

lasso_pipe.fit(X_train, y_train)
Pipeline(steps=[('preprocessor',
                 ColumnTransformer(transformers=[('num',
                                                  Pipeline(steps=[('imputer',
                                                                   SimpleImputer(strategy='median')),
                                                                  ('scaler',
                                                                   StandardScaler())]),
                                                  <sklearn.compose._column_transformer.make_column_selector object at 0x00000222F09495B0>),
                                                 ('cat',
                                                  Pipeline(steps=[('imputer',
                                                                   SimpleImputer(fill_value='missing',
                                                                                 strategy='constant')),
                                                                  ('onehot',
                                                                   OneHotEncoder(handle_unknown='ignore'))]),
                                                  <sklearn.compose._column_transformer.make_column_selector object at 0x00000222F0949940>)])),
                ('lasso', Lasso(alpha=1, max_iter=10000))])
# Metrics trainings set
pred_train = lasso_pipe.predict(X_train)
print('R squared training set:', round(skm.r2_score(y_train, pred_train)*100, 2))
print('MSE training set:', round(skm.mean_squared_error(y_train, pred_train), 2))
print('RMSE training set:', round(skm.mean_squared_error(y_train, pred_train, squared=False), 2))
R squared training set: 74.47
MSE training set: 2777618427.33
RMSE training set: 52703.12

R-squared der Lasso Regression hat sich im Vergleich zur Simple Regression im Modell zuvor nicht verändert. Auch MSE und RMSE sind gleich geblieben. Dies ist auf den initial gewählten, sehr kleinen Wert für alpha von 1 zurückzuführen. Ein geringer Wert in Alpha verändert die Höhe der Koeffizienten im Modell kaum. Daher bleiben auch die Berechnungen des Modells unverändert und somit die Metriken.

# Metrics validation set
pred = lasso_pipe.predict(X_test)
print('R squared validation set:', round(skm.r2_score(y_test, pred)*100, 2))
print('MSE validation set:', round(skm.mean_squared_error(y_test, pred), 2))
print('RMSE validation set:', round(skm.mean_squared_error(y_test, pred, squared=False), 2))
R squared validation set: 74.67
MSE validation set: 2877369775.47
RMSE validation set: 53641.12

Auch die Metriken des Validation Set haben sich nicht verändert im Vergleich zur Simple Regression.

# Obtain model coefficients and names
print('Koeffizienten und Namen: \n', list(zip(lasso_pipe.named_steps['lasso'].coef_, feature_names)))
Koeffizienten und Namen: 
 [(84852.14248638194, 'median_income'), (18734.010558726382, 'bedrooms_per_rooms'), (-19962.92053743495, 'people_per_household'), (11487.952227452783, 'housing_median_age'), (-57171.981131734494, 'x0_INLAND'), (3.3287087880073534e-11, 'x0_NEAR COAST')]
# get absolute values of coefficients
importance = np.abs(lasso_pipe.named_steps['lasso'].coef_)

sns.barplot(x=importance, 
            y=feature_names);
_images/regression_200_0.png

Das Diagramm visualiert die Feature Importance. Je länger der Balken, desto eine höhere Gewichtung bekommt das Feature im Modell. Die Feature Importance ergibt sich auch den Koeffizienten der Features. median_income ist das wichtigste Feature im Modell, wie an dem großen Balken sowie dem Koeffizienten des Features zu erkennen.
Die Ausprägung NEAR COAST der Variablen ocean_proximity wurde von der Lasso Regression aus der Regression fast entfernt. Dieses Verhalten ist bei Lasso gewünscht. Durch die L1-Regularisierung mit alpha führt der Algorithmus automatisch eine Feature Selection durch. Dabei werden relevante Features durch höhere Koeffizienten stärker gewichtet und unwichtige Features durch kleinere Koeffizienten schwächer gewichtet bzw. sogar komplett aus dem Modell entfernt, indem der Koeffizient auf den Wert 0 reduziert wird. Durch diese Methode soll ein overfitting des Modells durch zu viele Variablen verhindert werden.

k-fold cross validation zur Ermittlung von alpha

from sklearn.linear_model import LassoCV

# Lasso with 5 fold cross-validation
model_pipe = Pipeline(steps=[
    ('preprocessor', preprocessor),
    ('lassoCV', LassoCV(cv=5, random_state=0, max_iter=10000))
                        ])
# Fit model
model_pipe.fit(X_train, y_train)
Pipeline(steps=[('preprocessor',
                 ColumnTransformer(transformers=[('num',
                                                  Pipeline(steps=[('imputer',
                                                                   SimpleImputer(strategy='median')),
                                                                  ('scaler',
                                                                   StandardScaler())]),
                                                  <sklearn.compose._column_transformer.make_column_selector object at 0x00000222F09495B0>),
                                                 ('cat',
                                                  Pipeline(steps=[('imputer',
                                                                   SimpleImputer(fill_value='missing',
                                                                                 strategy='constant')),
                                                                  ('onehot',
                                                                   OneHotEncoder(handle_unknown='ignore'))]),
                                                  <sklearn.compose._column_transformer.make_column_selector object at 0x00000222F0949940>)])),
                ('lassoCV', LassoCV(cv=5, max_iter=10000, random_state=0))])
#Show best value of penalization chosen by cross validation:
alpha = model_pipe.named_steps['lassoCV'].alpha_
print('Best alpha chosen by cross validation:', alpha)
Best alpha chosen by cross validation: 79.0619710978722

Auf Basis von k-fold cross validation wurde der beste Wert für Alpha ermittelt. Dieser liegt bei 79,06. Bei k-fold cross validation wird das Trainingsdatenset nochmal in k kleinere Datensets unterteilt. Eines dieser kleineren Datensets wird zur Seite gelegt und das Modell mithilfe der anderen Datensätzen und einem Wert für Alpha trainiert. Anschließend wird mit dem zur Seite gelegten kleinen Datenset eine Validierung des Modells mit den unbekannten Daten durchgeführt.
Dieser Vorgang wird für jedes der k - Datensets für verschiedene Werte von Alpha wiederholt. Das Alpha mit dem besten Ergebnis bei der Validierung wird hier ausgegeben.

Mit dem zuvor ermittelten, besten Wert von Alpha wird nun noch ein finales Lasso Modell auf Basis aller Trainingsdaten trainiert.

# Set best alpha
lasso_best_pipe = Pipeline(steps=[
    ('preprocessor', preprocessor),
    ('lasso_best', Lasso(alpha=alpha))
                        ])
lasso_best_pipe.fit(X_train, y_train)
lasso_best = lasso_best_pipe.named_steps['lasso_best']
print('Model coefficients and names: \n', list(zip(lasso_best.coef_, feature_names)))
Model coefficients and names: 
 [(84682.01192725677, 'median_income'), (18560.71403670962, 'bedrooms_per_rooms'), (-19890.4973463913, 'people_per_household'), (11429.85201175577, 'housing_median_age'), (-56979.210567451206, 'x0_INLAND'), (5.031471360334191e-11, 'x0_NEAR COAST')]

Durch den größeren Wert von Alpha wurden die Koeffizienten der einzelnen Features etwas verkleinert. median_income hatte zuvor bei alpha=1 beispielsweise einen Wert von 84852,14. Nun liegt der Wert bei 84682,01.

# get absolute values of coefficients
importance = np.abs(lasso_best.coef_)

sns.barplot(x=importance, 
            y=feature_names);
_images/regression_210_0.png

An der Reihenfolge der Feature Importance hat sich durch die Veränderung von Alpha wenig verändert.

# Metrics trainings set
pred_train = lasso_best_pipe.predict(X_train)
print('R squared training set:', round(skm.r2_score(y_train, pred_train)*100, 2))
print('MSE training set:', round(skm.mean_squared_error(y_train, pred_train), 2))
print('RMSE training set:', round(skm.mean_squared_error(y_train, pred_train, squared=False), 2))
R squared training set: 74.47
MSE training set: 2777671055.36
RMSE training set: 52703.62
# Metrics tests set
pred = lasso_best_pipe.predict(X_test)
print('R squared validation set:', round(skm.r2_score(y_test, pred)*100, 2))
print('MSE validation set:', round(skm.mean_squared_error(y_test, pred), 2))
print('RMSE validation set:', round(skm.mean_squared_error(y_test, pred, squared=False), 2))
R squared validation set: 75.03
MSE validation set: 2817208350.85
RMSE validation set: 53077.38

Da sich die Koeffizienten durch den neuen Wert von Alpha kaum verändert haben, sind auch die Kennzahlen zur Performance Messung des Modells unverändert.

Model 7 - Natural Splines

In diesem Modell wird eine scikit learn spline in Kombination mit einer Data Preprocessing Pipeline trainiert.

# use a spline wit 4 knots and 3 degrees with a ridge regressions
spline = make_pipeline(SplineTransformer(n_knots=4, degree=3), 
                       Ridge(alpha=1))

# Integrate Preprocessor Pipeline
spline_pipe = Pipeline(steps=[
    ('preprocessor', preprocessor),
    ('spline', spline)
                        ])
                     
spline_pipe.fit(X_train, y_train)
Pipeline(steps=[('preprocessor',
                 ColumnTransformer(transformers=[('num',
                                                  Pipeline(steps=[('imputer',
                                                                   SimpleImputer(strategy='median')),
                                                                  ('scaler',
                                                                   StandardScaler())]),
                                                  <sklearn.compose._column_transformer.make_column_selector object at 0x00000222F09495B0>),
                                                 ('cat',
                                                  Pipeline(steps=[('imputer',
                                                                   SimpleImputer(fill_value='missing',
                                                                                 strategy='constant')),
                                                                  ('onehot',
                                                                   OneHotEncoder(handle_unknown='ignore'))]),
                                                  <sklearn.compose._column_transformer.make_column_selector object at 0x00000222F0949940>)])),
                ('spline',
                 Pipeline(steps=[('splinetransformer',
                                  SplineTransformer(n_knots=4)),
                                 ('ridge', Ridge(alpha=1))]))])
# Metrics trainings set
pred_train = spline_pipe.predict(X_train)
print('R squared training set:', round(skm.r2_score(y_train, pred_train)*100, 2))
print('MSE training set:', round(skm.mean_squared_error(y_train, pred_train), 2))
print('RMSE training set:', round(skm.mean_squared_error(y_train, pred_train, squared=False), 2))
R squared training set: 75.36
MSE training set: 2680356009.84
RMSE training set: 51772.15
# Metrics tests set
pred = spline_pipe.predict(X_test)
print('R squared validation set:', round(skm.r2_score(y_test, pred)*100, 2))
print('MSE validation set:', round(skm.mean_squared_error(y_test, pred), 2))
print('RMSE validation set:', round(skm.mean_squared_error(y_test, pred, squared=False), 2))
R squared validation set: 76.19
MSE validation set: 2687000792.6
RMSE validation set: 51836.29

Durch die Verwendung der Spline hat sich das Modell in Bezug auf R-sqared und auch RMSE nochmal verbessert. Das Modell deckt nun 75,36% der Varianz (Testdaten) bzw. 76,19% der Validationsdaten ab. Die durchschnittliche Abweichung der tatsächlichen Werten zu den berechneten Werten beträgt ca. 51800$.

# Create observations
xp = pd.DataFrame(data=np.linspace(X_test.min(),X_test.max(), 100), columns = X_train.drop(['ocean_proximity'], axis=1).columns)
xp = xp.assign(ocean_proximity=lambda xp: 'INLAND')


# Make some predictions
pred_new = spline_pipe.predict(xp)

# plot
sns.scatterplot(x=X_train['median_income'], y=y_train)

plt.plot(xp['median_income'], pred_new, label='Cubic spline with degree=3', color='orange')
plt.legend();
C:\Users\Admin\AppData\Local\Temp/ipykernel_5788/3121680389.py:2: FutureWarning:

Dropping of nuisance columns in DataFrame reductions (with 'numeric_only=None') is deprecated; in a future version this will raise TypeError.  Select only valid columns before calling the reduction.
_images/regression_221_1.png

Evaluierung der besten 2 Modelle mit Test Datenset

In diesem Kapitel wird eine Evaluierung der zwei besten Modellen mithilfe der zu Beginn beiseite gelegten Testdaten durchgeführt.
Zunächst müssen die Testdaten transformiert werden, sodass sie in gleicher Form wie die Trainingsdaten vorliegen und von den Modellen verwendet werden können.

# write out a dict with the mapping of old to new
remap_cat_dict = {
    'NEAR OCEAN': 'NEAR COAST',
    'NEAR BAY': 'NEAR COAST',
    'INLAND': 'INLAND',
    '<1H OCEAN': 'NEAR COAST',
    'ISLAND' : 'NEAR COAST'}

test_dataset.ocean_proximity = test_dataset.ocean_proximity.map(remap_cat_dict).astype('category')
y_test = test_dataset['median_house_value']

X_test = test_dataset[['median_income', 'ocean_proximity', 'bedrooms_per_rooms', 'people_per_household', 'housing_median_age']]

Evaluation Statsmodels OLS lm2

# Metrics tests set
X_test = X_test.assign(pred_test1=lambda X_test: lm2.predict(X_test))
r2_test1 = round(skm.r2_score(y_test, X_test.pred_test1)*100, 2)
print('R squared Statsmodels OLS with Test Set:', r2_test1)
print('MSE Statsmodels OLS with Test Set:', round(skm.mean_squared_error(y_test, X_test.pred_test1), 2))
print('RMSE Statsmodels OLS with Test Set:', round(skm.mean_squared_error(y_test, X_test.pred_test1, squared=False), 2))
R squared Statsmodels OLS with Test Set: -10.35
MSE Statsmodels OLS with Test Set: 14459492747.15
RMSE Statsmodels OLS with Test Set: 120247.63

R-squared ist in Model lm2 negativ. Somit generalisiert das Modell schlecht auf unbekannte Testdaten und stellt sogar eine Verschlechterung im Vergleich zur Verwendung des Durchschnitts als Modell dar.
Grund dafür könnte ein Overfitting während dem Modellierungsprozess sein.
Root Mean Squared Error und somit eine durchschnittliche Abweichung von tatsächlichem und berechnetem median_house_value in Höhe von 120247,63,95$ ist ebenso kein zufriedenstellendes Ergebnis.

Die berechneten Werte sollen in einem Scatterplot den tatsächlichen Werten gegenübergestellt werden. Aufgrund der Datenmenge sollen nur einige Datensätze visualisiert werden. Mithilfe eines Data Splittings werden zufällig einige Datenpunkte ausgewählt.

X_notplot, X_plot, y_notplot, y_plot = train_test_split(X_test, y_test, test_size=0.1, random_state=10)
# plot
sns.scatterplot(x=X_plot['median_income'], y=y_plot, label='actual')

plt.plot(X_plot['median_income'], X_plot['pred_test1'], label='predicted', color='orange', marker='.', linestyle='')
plt.legend();
_images/regression_232_0.png

Im Plot ist zu erkennen, dass bei bei einem median_income von einen Datenpunkt gibt, dessen median_house_value negativ geschätzt wurde. Dieser Wert ist unplausibel, da ein Haus keinen negativen Wert hat.
Den Koeffizienten des Modells lm2 ist zu entnehmen, dass nur das Feature people_per_household eine negative Abhängigkeit auf median_house_value hat. Die Ursache für einen negativen Wert wird also durch einen hohen Wert im Feature people_per_household verursacht.

lm2.params
Intercept                        -49511.050850
ocean_proximity[T.NEAR COAST]     57629.618091
median_income                     48571.655721
bedrooms_per_rooms               363707.323412
people_per_household             -27523.527384
housing_median_age                  923.267802
dtype: float64
test_dataset.loc[test_dataset['people_per_household'] > 100]
longitude latitude housing_median_age total_rooms total_bedrooms population households median_income median_house_value ocean_proximity people_per_household rooms_per_household bedrooms_per_people bedrooms_per_rooms
13034 -121.15 38.69 52 240 44 6675 29 6.1359 225000.0 INLAND 230.172414 8.275862 0.006592 0.183333

Tatsächlich existiert ein Datensatz mit über 100 people_per_household. Dieser ist vermutlich für den negativen R-squared Wert verantwortlich. Daher wird dieser testweise entfernt und R-squared erneut berechnet.

X_test2 = X_test.drop(13034)
y_test2 = y_test.drop(13034)
# Metrics tests set
X_test2 = X_test2.assign(pred_test1=lambda X_test2: lm2.predict(X_test2))
r2_test1 = round(skm.r2_score(y_test2, X_test2.pred_test1)*100, 2)
print('R squared Statsmodels OLS with Test Set:', r2_test1)
print('MSE Statsmodels OLS with Test Set:', round(skm.mean_squared_error(y_test2, X_test2.pred_test1), 2))
print('RMSE Statsmodels OLS with Test Set:', round(skm.mean_squared_error(y_test2, X_test2.pred_test1, squared=False), 2))
R squared Statsmodels OLS with Test Set: 61.36
MSE Statsmodels OLS with Test Set: 5064559536.22
RMSE Statsmodels OLS with Test Set: 71165.72

Wie vermutet, lag der Grund für den negativen Wert in diesem Ausreiser. Jedoch darf auf den Testdaten gemäß der Literatur keine EDA oder Feature Engineering betrieben werden.
Das Modell lm2 ist daher kein gutes Modell, welches ausreichend unbekannte Daten generalisiert.

Evaluation Scikit Learn Spline

#Calculate R2
X_test = X_test.assign(pred_test2=lambda X_test: spline_pipe.predict(X_test))
r2_test = round(skm.r2_score(y_test, X_test.pred_test2)*100, 2)
print('R squared test set:', r2_test)
mse_test =skm.mean_squared_error(y_test, X_test.pred_test2)
print('MSE test set', round(mse_test, 2))
print('RMSE Statsmodels OLS with Test Set:', round(skm.mean_squared_error(y_test, X_test.pred_test2, squared=False), 2))
R squared test set: 65.54
MSE test set 4515256968.94
RMSE Statsmodels OLS with Test Set: 67195.66

Dieses Modell berechnet aus dem kompletten Test Datenset einen R-squared Wert von 65,54%. Damit wird das zu Beginn gesetzte Ziel, einen R-squared von 70% zu erreichen unter Verwendung der Testdaten leider verfehlt. Unter Verwendung der Trainingsdaten wurde der Wert jedoch erreicht.
RMSE erzielt einen Wert von 67195,66$ und ist damit auch eher hoch.

X_plot = X_plot.assign(pred_test2=lambda X_test: spline_pipe.predict(X_test))
# plot
sns.scatterplot(x=X_plot['median_income'], y=y_plot, label='actual')

plt.plot(X_plot['median_income'], X_plot['pred_test2'], label='predicted', color='orange', marker='.', linestyle='')
plt.legend();
_images/regression_243_0.png

Den Analysen entsprechend ist das beste Modell ist die Spline von Scikit Learn mit einem R-squared der Testdaten von 65,25%.
Gemäß dem Data-Science Lifecycle müsste dieses Modell jetzt deployed werden. Dieser Schritt wird allerdings nicht mehr im Rahmen dieses Projekts durchgeführt.